source: sasmodels/sasmodels/kernel_iq.cl @ 8698a0d

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 8698a0d was 8698a0d, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

revise api for oriented shapes, allowing jitter in the frame of the sample
  • Property mode set to 100644
File size: 14.9 KB
Line 
1
2/*
3    ##########################################################
4    #                                                        #
5    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
6    #   !!                                              !!   #
7    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
8    #   !!                                              !!   #
9    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
10    #                                                        #
11    ##########################################################
12*/
13
14#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
15#define _PAR_BLOCK_
16
17typedef struct {
18#if MAX_PD > 0
19    int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable
20    int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector
21    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
22    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
23#endif // MAX_PD > 0
24    int32_t num_eval;           // total number of voxels in hypercube
25    int32_t num_weights;        // total length of the weights vector
26    int32_t num_active;         // number of non-trivial pd loops
27    int32_t theta_par;          // id of spherical correction variable (not used)
28} ProblemDetails;
29
30// Intel HD 4000 needs private arrays to be a multiple of 4 long
31typedef struct {
32    PARAMETER_TABLE
33} ParameterTable;
34typedef union {
35    ParameterTable table;
36    double vector[4*((NUM_PARS+3)/4)];
37} ParameterBlock;
38#endif // _PAR_BLOCK_
39
40
41#if defined(MAGNETIC) && NUM_MAGNETIC>0
42
43// Return value restricted between low and high
44static double clip(double value, double low, double high)
45{
46  return (value < low ? low : (value > high ? high : value));
47}
48
49// Compute spin cross sections given in_spin and out_spin
50// To convert spin cross sections to sld b:
51//     uu * (sld - m_sigma_x);
52//     dd * (sld + m_sigma_x);
53//     ud * (m_sigma_y + 1j*m_sigma_z);
54//     du * (m_sigma_y - 1j*m_sigma_z);
55static void set_spins(double in_spin, double out_spin, double spins[4])
56{
57  in_spin = clip(in_spin, 0.0, 1.0);
58  out_spin = clip(out_spin, 0.0, 1.0);
59  spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd
60  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du
61  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud
62  spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu
63}
64
65static double mag_sld(double qx, double qy, double p,
66                       double mx, double my, double sld)
67{
68    const double perp = qy*mx - qx*my;
69    return sld + perp*p;
70}
71
72//#endif // MAGNETIC
73
74// TODO: way too hackish
75// For the 1D kernel, CALL_IQ_[A,AC,ABC] and MAGNETIC are not defined
76// so view_direct *IS NOT* included
77// For the 2D kernel, CALL_IQ_[A,AC,ABC] is defined but MAGNETIC is not
78// so view_direct *IS* included
79// For the magnetic kernel, CALL_IQ_[A,AC,ABC] is defined, but so is MAGNETIC
80// so view_direct *IS NOT* included
81#else // !MAGNETIC
82
83// ===== Implement jitter in orientation =====
84// To change the definition of the angles, run explore/angles.py, which
85// uses sympy to generate the equations.
86
87#if defined(CALL_IQ_AC) // oriented symmetric
88static double
89view_direct(double qx, double qy,
90             double theta, double phi,
91             ParameterTable table)
92{
93    double sin_theta, cos_theta;
94    double sin_phi, cos_phi;
95
96    // reverse view
97    SINCOS(theta*M_PI_180, sin_theta, cos_theta);
98    SINCOS(phi*M_PI_180, sin_phi, cos_phi);
99    const double qa = qx*cos_phi*cos_theta + qy*sin_phi*cos_theta;
100    const double qb = -qx*sin_phi + qy*cos_phi;
101    const double qc = qx*sin_theta*cos_phi + qy*sin_phi*sin_theta;
102
103    // reverse jitter after view
104    SINCOS(table.theta*M_PI_180, sin_theta, cos_theta);
105    SINCOS(table.phi*M_PI_180, sin_phi, cos_phi);
106    const double dqc = qa*sin_theta - qb*sin_phi*cos_theta + qc*cos_phi*cos_theta;
107
108    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2
109    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy);
110
111    return CALL_IQ_AC(dqa, dqc, table);
112}
113
114#elif defined(CALL_IQ_ABC) // oriented asymmetric
115
116static double
117view_direct(double qx, double qy,
118             double theta, double phi, double psi,
119             ParameterTable table)
120{
121    double sin_theta, cos_theta;
122    double sin_phi, cos_phi;
123    double sin_psi, cos_psi;
124
125    // reverse view
126    SINCOS(theta*M_PI_180, sin_theta, cos_theta);
127    SINCOS(phi*M_PI_180, sin_phi, cos_phi);
128    SINCOS(psi*M_PI_180, sin_psi, cos_psi);
129    const double qa = qx*(sin_phi*sin_psi + cos_phi*cos_psi*cos_theta) + qy*(sin_phi*cos_psi*cos_theta - sin_psi*cos_phi);
130    const double qb = qx*(-sin_phi*cos_psi + sin_psi*cos_phi*cos_theta) + qy*(sin_phi*sin_psi*cos_theta + cos_phi*cos_psi);
131    const double qc = qx*sin_theta*cos_phi + qy*sin_phi*sin_theta;
132
133    // reverse jitter after view
134    SINCOS(table.theta*M_PI_180, sin_theta, cos_theta);
135    SINCOS(table.phi*M_PI_180, sin_phi, cos_phi);
136    SINCOS(table.psi*M_PI_180, sin_psi, cos_psi);
137    const double dqa = qa*cos_psi*cos_theta + qb*(sin_phi*sin_theta*cos_psi - sin_psi*cos_phi) + qc*(-sin_phi*sin_psi - sin_theta*cos_phi*cos_psi);
138    const double dqb = qa*sin_psi*cos_theta + qb*(sin_phi*sin_psi*sin_theta + cos_phi*cos_psi) + qc*(sin_phi*cos_psi - sin_psi*sin_theta*cos_phi);
139    const double dqc = qa*sin_theta - qb*sin_phi*cos_theta + qc*cos_phi*cos_theta;
140
141    return CALL_IQ_ABC(dqa, dqb, dqc, table);
142}
143#endif
144
145#endif // !MAGNETIC
146
147kernel
148void KERNEL_NAME(
149    int32_t nq,                 // number of q values
150    const int32_t pd_start,     // where we are in the polydispersity loop
151    const int32_t pd_stop,      // where we are stopping in the polydispersity loop
152    global const ProblemDetails *details,
153    global const double *values,
154    global const double *q, // nq q values, with padding to boundary
155    global double *result,  // nq+1 return values, again with padding
156    const double cutoff     // cutoff in the polydispersity weight product
157    )
158{
159
160  // who we are and what element we are working with
161  const int q_index = get_global_id(0);
162  if (q_index >= nq) return;
163
164  // Storage for the current parameter values.  These will be updated as we
165  // walk the polydispersity cube.
166  ParameterBlock local_values;
167
168#if defined(MAGNETIC) && NUM_MAGNETIC>0
169  // Location of the sld parameters in the parameter vector.
170  // These parameters are updated with the effective sld due to magnetism.
171  #if NUM_MAGNETIC > 3
172  const int32_t slds[] = { MAGNETIC_PARS };
173  #endif
174
175  // TODO: could precompute these outside of the kernel.
176  // Interpret polarization cross section.
177  //     up_frac_i = values[NUM_PARS+2];
178  //     up_frac_f = values[NUM_PARS+3];
179  //     up_angle = values[NUM_PARS+4];
180  double spins[4];
181  double cos_mspin, sin_mspin;
182  set_spins(values[NUM_PARS+2], values[NUM_PARS+3], spins);
183  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
184#endif // MAGNETIC
185
186#if defined(CALL_IQ_AC) // oriented symmetric
187  const double theta = values[details->theta_par+2];
188  const double phi = values[details->theta_par+3];
189#elif defined(CALL_IQ_ABC) // oriented asymmetric
190  const double theta = values[details->theta_par+2];
191  const double phi = values[details->theta_par+3];
192  const double psi = values[details->theta_par+4];
193#endif
194
195  // Fill in the initial variables
196  //   values[0] is scale
197  //   values[1] is background
198  for (int i=0; i < NUM_PARS; i++) {
199    local_values.vector[i] = values[2+i];
200//if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]);
201  }
202//if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
203//if (q_index==0) printf("start:%d stop:%d\n", pd_start, pd_stop);
204
205  double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
206  double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
207//if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, this_result);
208
209#if MAX_PD>0
210  global const double *pd_value = values + NUM_VALUES;
211  global const double *pd_weight = pd_value + details->num_weights;
212#endif
213
214  // Jump into the middle of the polydispersity loop
215#if MAX_PD>4
216  int n4=details->pd_length[4];
217  int i4=(pd_start/details->pd_stride[4])%n4;
218  const int p4=details->pd_par[4];
219  global const double *v4 = pd_value + details->pd_offset[4];
220  global const double *w4 = pd_weight + details->pd_offset[4];
221#endif
222#if MAX_PD>3
223  int n3=details->pd_length[3];
224  int i3=(pd_start/details->pd_stride[3])%n3;
225  const int p3=details->pd_par[3];
226  global const double *v3 = pd_value + details->pd_offset[3];
227  global const double *w3 = pd_weight + details->pd_offset[3];
228//if (q_index==0) printf("offset %d: %d %d\n", 3, details->pd_offset[3], NUM_VALUES);
229#endif
230#if MAX_PD>2
231  int n2=details->pd_length[2];
232  int i2=(pd_start/details->pd_stride[2])%n2;
233  const int p2=details->pd_par[2];
234  global const double *v2 = pd_value + details->pd_offset[2];
235  global const double *w2 = pd_weight + details->pd_offset[2];
236#endif
237#if MAX_PD>1
238  int n1=details->pd_length[1];
239  int i1=(pd_start/details->pd_stride[1])%n1;
240  const int p1=details->pd_par[1];
241  global const double *v1 = pd_value + details->pd_offset[1];
242  global const double *w1 = pd_weight + details->pd_offset[1];
243#endif
244#if MAX_PD>0
245  int n0=details->pd_length[0];
246  int i0=(pd_start/details->pd_stride[0])%n0;
247  const int p0=details->pd_par[0];
248  global const double *v0 = pd_value + details->pd_offset[0];
249  global const double *w0 = pd_weight + details->pd_offset[0];
250#endif
251
252
253  int step = pd_start;
254
255
256#if MAX_PD>4
257  const double weight5 = 1.0;
258  while (i4 < n4) {
259    local_values.vector[p4] = v4[i4];
260    double weight4 = w4[i4] * weight5;
261//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, local_values.vector[p4], weight4);
262#elif MAX_PD>3
263    const double weight4 = 1.0;
264#endif
265#if MAX_PD>3
266  while (i3 < n3) {
267    local_values.vector[p3] = v3[i3];
268    double weight3 = w3[i3] * weight4;
269//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, local_values.vector[p3], weight3);
270#elif MAX_PD>2
271    const double weight3 = 1.0;
272#endif
273#if MAX_PD>2
274  while (i2 < n2) {
275    local_values.vector[p2] = v2[i2];
276    double weight2 = w2[i2] * weight3;
277//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, local_values.vector[p2], weight2);
278#elif MAX_PD>1
279    const double weight2 = 1.0;
280#endif
281#if MAX_PD>1
282  while (i1 < n1) {
283    local_values.vector[p1] = v1[i1];
284    double weight1 = w1[i1] * weight2;
285//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, local_values.vector[p1], weight1);
286#elif MAX_PD>0
287    const double weight1 = 1.0;
288#endif
289#if MAX_PD>0
290  while(i0 < n0) {
291    local_values.vector[p0] = v0[i0];
292    double weight0 = w0[i0] * weight1;
293//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0);
294#else
295    const double weight0 = 1.0;
296#endif
297
298//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n"); }
299//#if defined(CALL_IQ_AC) // oriented symmetric
300//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); printf("theta=%g dtheta=%g",theta,local_values.table.theta); printf("\n"); }
301//#endif
302
303    #ifdef INVALID
304    if (!INVALID(local_values.table))
305    #endif
306    {
307      // Accumulate I(q)
308      // Note: weight==0 must always be excluded
309      if (weight0 > cutoff) {
310        pd_norm += weight0 * CALL_VOLUME(local_values.table);
311
312#if defined(MAGNETIC) && NUM_MAGNETIC > 0
313        const double qx = q[2*q_index];
314        const double qy = q[2*q_index+1];
315        const double qsq = qx*qx + qy*qy;
316
317        // Constant across orientation, polydispersity for given qx, qy
318        double scattering = 0.0;
319        // TODO: what is the magnetic scattering at q=0
320        if (qsq > 1.e-16) {
321          double p[4];  // dd, du, ud, uu
322          p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq;
323          p[3] = -p[0];
324          p[1] = p[2] = (qy*sin_mspin - qx*cos_mspin)/qsq;
325
326          for (int index=0; index<4; index++) {
327            const double xs = spins[index];
328            if (xs > 1.e-8) {
329              const int spin_flip = (index==1) || (index==2);
330              const double pk = p[index];
331              for (int axis=0; axis<=spin_flip; axis++) {
332                #define M1 NUM_PARS+5
333                #define M2 NUM_PARS+8
334                #define M3 NUM_PARS+13
335                #define SLD(_M_offset, _sld_offset) \
336                    local_values.vector[_sld_offset] = xs * (axis \
337                    ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \
338                    : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \
339                              (spin_flip ? 0.0 : values[_sld_offset+2])))
340                #if NUM_MAGNETIC==1
341                    SLD(M1, MAGNETIC_PAR1);
342                #elif NUM_MAGNETIC==2
343                    SLD(M1, MAGNETIC_PAR1);
344                    SLD(M2, MAGNETIC_PAR2);
345                #elif NUM_MAGNETIC==3
346                    SLD(M1, MAGNETIC_PAR1);
347                    SLD(M2, MAGNETIC_PAR2);
348                    SLD(M3, MAGNETIC_PAR3);
349                #else
350                for (int sk=0; sk<NUM_MAGNETIC; sk++) {
351                    SLD(M1+3*sk, slds[sk]);
352                }
353                #endif
354#  if defined(CALL_IQ_A) // unoriented
355                scattering += CALL_IQ_A(sqrt(qsq), local_values.table);
356#  elif defined(CALL_IQ_AC) // oriented symmetric
357                scattering += view_direct(qx, qy, theta, phi, local_values.table);
358#  elif defined(CALL_IQ_ABC) // oriented asymmetric
359                scattering += view_direct(qx, qy, theta, phi, psi, local_values.table);
360#  endif
361              }
362            }
363          }
364        }
365#elif defined(CALL_IQ) // 1d, not MAGNETIC
366          const double scattering = CALL_IQ(q[q_index], local_values.table);
367#else  // 2d data, not magnetic
368          const double qx = q[2*q_index];
369          const double qy = q[2*q_index+1];
370#  if defined(CALL_IQ_A) // unoriented
371          const double scattering = CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table);
372#  elif defined(CALL_IQ_AC) // oriented symmetric
373          const double scattering = view_direct(qx, qy, theta, phi, local_values.table);
374#  elif defined(CALL_IQ_ABC) // oriented asymmetric
375          const double scattering = view_direct(qx, qy, theta, phi, psi, local_values.table);
376#  endif
377#endif // !MAGNETIC
378        this_result += weight0 * scattering;
379      }
380    }
381    ++step;
382#if MAX_PD>0
383    if (step >= pd_stop) break;
384    ++i0;
385  }
386  i0 = 0;
387#endif
388#if MAX_PD>1
389    if (step >= pd_stop) break;
390    ++i1;
391  }
392  i1 = 0;
393#endif
394#if MAX_PD>2
395    if (step >= pd_stop) break;
396    ++i2;
397  }
398  i2 = 0;
399#endif
400#if MAX_PD>3
401    if (step >= pd_stop) break;
402    ++i3;
403  }
404  i3 = 0;
405#endif
406#if MAX_PD>4
407    if (step >= pd_stop) break;
408    ++i4;
409  }
410  i4 = 0;
411#endif
412
413//if (q_index==0) printf("res: %g/%g\n", this_result, pd_norm);
414  // Remember the current result and the updated norm.
415  result[q_index] = this_result;
416  if (q_index == 0) result[nq] = pd_norm;
417}
Note: See TracBrowser for help on using the repository browser.