Changes in / [e9b17b18:7e923c2] in sasmodels


Ignore:
Files:
45 added
13 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/core.py

    r2dcd6e7 rc036ddb  
    205205 
    206206    numpy_dtype, fast, platform = parse_dtype(model_info, dtype, platform) 
    207  
    208207    source = generate.make_source(model_info) 
    209208    if platform == "dll": 
    … …  
    265264    # Assign default platform, overriding ocl with dll if OpenCL is unavailable 
    266265    # If opencl=False OpenCL is switched off 
    267  
    268266    if platform is None: 
    269267        platform = "ocl" 

  • sasmodels/data.py

    r581661f rc036ddb  
    329329    else: 
    330330        dq = resolution * q 
    331  
    332331    data = Data1D(q, Iq, dx=dq, dy=dIq) 
    333332    data.filename = "fake data" 

  • sasmodels/generate.py

    rd86f0fc rc036ddb  
    671671 
    672672 
    673 # type in IQXY pattern could be single, float, double, long double, ... 
    674 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
     673_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]", 
    675674                           flags=re.MULTILINE) 
    676675def find_xy_mode(source): 
    … …  
    701700    return 'qa' 
    702701 
     702# Note: not presently used.  Need to know whether Fq is available before 
     703# trying to compile the source.  Leave the code here in case we decide that 
     704# define have_Fq for each form factor is too tedious and error prone. 
     705_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE) 
     706def contains_Fq(source): 
     707    # type: (List[str]) -> bool 
     708    """ 
     709    Return True if C source defines "void Fq(". 
     710    """ 
     711    for code in source: 
     712        m = _FQ_PATTERN.search(code) 
     713        if m is not None: 
     714            return True 
     715    return False 
    703716 
    704717def _add_source(source, code, path, lineno=1): 
    … …  
    730743    # dispersion.  Need to be careful that necessary parameters are available 
    731744    # for computing volume even if we allow non-disperse volume parameters. 
    732  
    733745    partable = model_info.parameters 
    734746 
    … …  
    743755    for path, code in user_code: 
    744756        _add_source(source, code, path) 
    745  
    746757    if model_info.c_code: 
    747758        _add_source(source, model_info.c_code, model_info.filename, 
    … …  
    751762    q, qx, qy, qab, qa, qb, qc \ 
    752763        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
     764 
    753765    # Generate form_volume function, etc. from body only 
    754766    if isinstance(model_info.form_volume, str): 
    … …  
    789801    source.append("\\\n".join(p.as_definition() 
    790802                              for p in partable.kernel_parameters)) 
    791  
    792803    # Define the function calls 
    793804    if partable.form_volume_parameters: 
    … …  
    800811        call_volume = "#define CALL_VOLUME(v) 1.0" 
    801812    source.append(call_volume) 
    802  
    803813    model_refs = _call_pars("_v.", partable.iq_parameters) 
    804     pars = ",".join(["_q"] + model_refs) 
    805     call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     814 
     815    if model_info.have_Fq: 
     816        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs) 
     817        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars 
     818        clear_iq = "#undef CALL_FQ" 
     819    else: 
     820        pars = ",".join(["_q"] + model_refs) 
     821        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     822        clear_iq = "#undef CALL_IQ" 
    806823    if xy_mode == 'qabc': 
    807824        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    … …  
    812829        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    813830        clear_iqxy = "#undef CALL_IQ_AC" 
    814     elif xy_mode == 'qa': 
     831    elif xy_mode == 'qa' and not model_info.have_Fq: 
    815832        pars = ",".join(["_qa"] + model_refs) 
    816833        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    817834        clear_iqxy = "#undef CALL_IQ_A" 
     835    elif xy_mode == 'qa' and model_info.have_Fq: 
     836        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs) 
     837        # Note: uses rare C construction (expr1, expr2) which computes 
     838        # expr1 then expr2 and evaluates to expr2.  This allows us to 
     839        # leave it looking like a function even though it is returning 
     840        # its values by reference. 
     841        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars 
     842        clear_iqxy = "#undef CALL_FQ_A" 
    818843    elif xy_mode == 'qxy': 
    819844        orientation_refs = _call_pars("_v.", partable.orientation_parameters) 
    … …  
    831856    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    832857               if p.type == 'sld'] 
    833  
    834858    # Fill in definitions for numbers of parameters 
    835859    source.append("#define MAX_PD %s"%partable.max_pd) 
    … …  
    839863    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    840864    source.append("#define PROJECTION %d"%PROJECTION) 
    841  
    842865    # TODO: allow mixed python/opencl kernels? 
    843  
    844     ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    845     dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
     866    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
     867    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
     868 
    846869    result = { 
    847870        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    848871        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    849872    } 
    850  
    851873    return result 
    852874 
    853875 
    854 def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
     876def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name): 
    855877    # type: ([str,str], str, str, str) -> List[str] 
    856878    code = kernel[0] 
    … …  
    862884        '#line 1 "%s Iq"' % path, 
    863885        code, 
    864         "#undef CALL_IQ", 
     886        clear_iq, 
    865887        "#undef KERNEL_NAME", 
    866888        ] 

  • sasmodels/kernel_iq.c

    r7c35fda rc036ddb  
    2929//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
    3030//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
     31//  CALL_FQ(q, F1, F2, table) : call the Fq function for 1D calcs. 
     32//  CALL_FQ_A(q, F1, F2, table) : call the Iq function with |q| for 2D data. 
    3133//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 
    3234//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes 
    … …  
    3638//  PROJECTION : equirectangular=1, sinusoidal=2 
    3739//      see explore/jitter.py for definitions. 
     40 
    3841 
    3942#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
    … …  
    8487  out_spin = clip(out_spin, 0.0, 1.0); 
    8588  // Previous version of this function took the square root of the weights, 
    86   // under the assumption that  
     89  // under the assumption that 
    8790  // 
    8891  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
    … …  
    270273#endif // _QABC_SECTION 
    271274 
    272  
    273275// ==================== KERNEL CODE ======================== 
    274  
     276#define COMPUTE_F1_F2 defined(CALL_FQ) 
    275277kernel 
    276278void KERNEL_NAME( 
    … …  
    341343  // version must loop over all q. 
    342344  #ifdef USE_OPENCL 
    343     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    344     double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
     345    #if COMPUTE_F1_F2 
     346      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     347      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     348      double this_F2 = (pd_start == 0 ? 0.0 : result[q_index]); 
     349      double this_F1 = (pd_start == 0 ? 0.0 : result[q_index+nq+1]); 
     350    #else 
     351      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     352      double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
     353    #endif 
    345354  #else // !USE_OPENCL 
    346     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     355    #if COMPUTE_F1_F2 
     356      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     357      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     358    #else 
     359      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     360    #endif 
    347361    if (pd_start == 0) { 
    348362      #ifdef USE_OPENMP 
    349363      #pragma omp parallel for 
    350364      #endif 
    351       for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     365      #if COMPUTE_F1_F2 
     366          for (int q_index=0; q_index < 2*nq+2; q_index++) result[q_index] = 0.0; 
     367      #else 
     368          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     369      #endif 
    352370    } 
    353371    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    … …  
    442460// inner loop and defines the macros that use them. 
    443461 
    444 #if defined(CALL_IQ) 
    445   // unoriented 1D 
     462 
     463#if defined(CALL_FQ) // COMPUTE_F1_F2 is true 
     464  // unoriented 1D returning <F> and <F^2> 
     465  double qk; 
     466  double F1, F2; 
     467  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     468  #define BUILD_ROTATION() do {} while(0) 
     469  #define APPLY_ROTATION() do {} while(0) 
     470  #define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table) 
     471 
     472#elif defined(CALL_FQ_A) 
     473  // unoriented 2D return <F> and <F^2> 
     474  double qx, qy; 
     475  double F1, F2; 
     476  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     477  #define BUILD_ROTATION() do {} while(0) 
     478  #define APPLY_ROTATION() do {} while(0) 
     479  #define CALL_KERNEL() CALL_FQ_A(sqrt(qx*qx+qy*qy),F1,F2,local_values.table) 
     480 
     481#elif defined(CALL_IQ) 
     482  // unoriented 1D return <F^2> 
    446483  double qk; 
    447484  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
    448485  #define BUILD_ROTATION() do {} while(0) 
    449486  #define APPLY_ROTATION() do {} while(0) 
    450   #define CALL_KERNEL() CALL_IQ(qk, local_values.table) 
     487  #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 
    451488 
    452489#elif defined(CALL_IQ_A) 
    … …  
    482519  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
    483520  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     521 
    484522#elif defined(CALL_IQ_XY) 
    485523  // direct call to qx,qy calculator 
    … …  
    495533// logic in the IQ_AC and IQ_ABC branches.  This will become more important 
    496534// if we implement more projections, or more complicated projections. 
    497 #if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation 
     535#if defined(CALL_IQ) || defined(CALL_IQ_A) || defined(CALL_FQ) || defined(CALL_FQ_A) 
     536  // no orientation 
    498537  #define APPLY_PROJECTION() const double weight=weight0 
    499538#elif defined(CALL_IQ_XY) // pass orientation to the model 
    … …  
    646685    if (weight > cutoff) { 
    647686      pd_norm += weight * CALL_VOLUME(local_values.table); 
     687      #if COMPUTE_F1_F2 
     688      weight_norm += weight; 
     689      #endif 
    648690      BUILD_ROTATION(); 
    649691 
    … …  
    693735          } 
    694736        #else  // !MAGNETIC 
    695           const double scattering = CALL_KERNEL(); 
     737          #if COMPUTE_F1_F2 
     738            CALL_KERNEL(); // sets F1 and F2 by reference 
     739          #else 
     740            const double scattering = CALL_KERNEL(); 
     741          #endif 
    696742        #endif // !MAGNETIC 
    697743//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 
    698744 
    699745        #ifdef USE_OPENCL 
    700           this_result += weight * scattering; 
     746          #if COMPUTE_F1_F2 
     747            this_F1 += weight * F1; 
     748            this_F2 += weight * F2; 
     749          #else 
     750            this_result += weight * scattering; 
     751          #endif 
    701752        #else // !USE_OPENCL 
    702           result[q_index] += weight * scattering; 
     753          #if COMPUTE_F1_F2 
     754            result[q_index] += weight * F2; 
     755            result[q_index+nq+1] += weight * F1; 
     756          #else 
     757            result[q_index] += weight * scattering; 
     758          #endif 
    703759        #endif // !USE_OPENCL 
    704760      } 
    705761    } 
    706762  } 
    707  
    708763// close nested loops 
    709764++step; 
    … …  
    726781// Remember the current result and the updated norm. 
    727782#ifdef USE_OPENCL 
    728   result[q_index] = this_result; 
    729   if (q_index == 0) result[nq] = pd_norm; 
     783  #if COMPUTE_F1_F2 
     784    result[q_index] = this_F2; 
     785    result[q_index+nq+1] = this_F1; 
     786    if (q_index == 0) { 
     787      result[nq] = pd_norm; 
     788      result[2*nq+1] = weight_norm; 
     789    } 
     790  #else 
     791    result[q_index] = this_result; 
     792    if (q_index == 0) result[nq] = pd_norm; 
     793  #endif 
     794 
    730795//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
    731796#else // !USE_OPENCL 
    732   result[nq] = pd_norm; 
     797  #if COMPUTE_F1_F2 
     798    result[nq] = pd_norm; 
     799    result[2*nq+1] = weight_norm; 
     800  #else 
     801    result[nq] = pd_norm; 
     802  #endif 
    733803//printf("res: %g/%g\n", result[0], pd_norm); 
    734804#endif // !USE_OPENCL 
    735805 
    736806// ** clear the macros in preparation for the next kernel ** 
     807#undef COMPUTE_F1_F2 
    737808#undef PD_INIT 
    738809#undef PD_OPEN 

  • sasmodels/kernelcl.py

    rd86f0fc rc036ddb  
    538538        self.dtype = dtype 
    539539        self.dim = '2d' if q_input.is_2d else '1d' 
    540         # plus three for the normalization values 
    541         self.result = np.empty(q_input.nq+1, dtype) 
     540        # leave room for f1/f2 results in case we need to compute beta for 1d models 
     541        num_returns = 1 if self.dim == '2d' else 2  # 
     542        # plus 1 for the normalization value 
     543        self.result = np.empty((q_input.nq+1)*num_returns, dtype) 
    542544 
    543545        # Inputs and outputs for each kernel call 
    … …  
    547549 
    548550        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    549                                   q_input.global_size[0] * dtype.itemsize) 
     551                                  q_input.global_size[0] * num_returns * dtype.itemsize) 
    550552        self.q_input = q_input # allocated by GpuInput above 
    551553 
    … …  
    556558                     else np.float32)  # will never get here, so use np.float32 
    557559 
    558     def __call__(self, call_details, values, cutoff, magnetic): 
     560    def Iq(self, call_details, values, cutoff, magnetic): 
     561        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
     562        self._call_kernel(call_details, values, cutoff, magnetic) 
     563        #print("returned",self.q_input.q, self.result) 
     564        pd_norm = self.result[self.q_input.nq] 
     565        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
     566        background = values[1] 
     567        #print("scale",scale,background) 
     568        return scale*self.result[:self.q_input.nq] + background 
     569    __call__ = Iq 
     570 
     571    def beta(self, call_details, values, cutoff, magnetic): 
     572        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
     573        if self.dim == '2d': 
     574            raise NotImplementedError("beta not yet supported for 2D") 
     575        self._call_kernel(call_details, values, cutoff, magnetic) 
     576        w_norm = self.result[2*self.q_input.nq + 1] 
     577        pd_norm = self.result[self.q_input.nq] 
     578        if w_norm == 0.: 
     579            w_norm = 1. 
     580        F2 = self.result[:self.q_input.nq]/w_norm 
     581        F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm 
     582        volume_avg = pd_norm/w_norm 
     583        return F1, F2, volume_avg 
     584 
     585    def _call_kernel(self, call_details, values, cutoff, magnetic): 
    559586        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    560587        context = self.queue.context 
    … …  
    573600        #print("Calling OpenCL") 
    574601        #call_details.show(values) 
    575         # Call kernel and retrieve results 
     602        #Call kernel and retrieve results 
    576603        wait_for = None 
    577604        last_nap = time.clock() 
    … …  
    598625                v.release() 
    599626 
    600         pd_norm = self.result[self.q_input.nq] 
    601         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    602         background = values[1] 
    603         #print("scale",scale,values[0],self.result[self.q_input.nq],background) 
    604         return scale*self.result[:self.q_input.nq] + background 
    605         # return self.result[:self.q_input.nq] 
    606  
    607627    def release(self): 
    608628        # type: () -> None 

  • sasmodels/kerneldll.py

    r1a3559f rc036ddb  
    375375    def __init__(self, kernel, model_info, q_input): 
    376376        # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 
     377        #,model_info,q_input) 
    377378        self.kernel = kernel 
    378379        self.info = model_info 
    … …  
    380381        self.dtype = q_input.dtype 
    381382        self.dim = '2d' if q_input.is_2d else '1d' 
    382         self.result = np.empty(q_input.nq+1, q_input.dtype) 
     383        # leave room for f1/f2 results in case we need to compute beta for 1d models 
     384        num_returns = 1 if self.dim == '2d' else 2  # 
     385        # plus 1 for the normalization value 
     386        self.result = np.empty((q_input.nq+1)*num_returns, self.dtype) 
    383387        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    384388                     else np.float64 if self.q_input.dtype == generate.F64 
    385389                     else np.float128) 
    386390 
    387     def __call__(self, call_details, values, cutoff, magnetic): 
     391    def Iq(self, call_details, values, cutoff, magnetic): 
    388392        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    389  
     393        self._call_kernel(call_details, values, cutoff, magnetic) 
     394        #print("returned",self.q_input.q, self.result) 
     395        pd_norm = self.result[self.q_input.nq] 
     396        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
     397        background = values[1] 
     398        #print("scale",scale,background) 
     399        return scale*self.result[:self.q_input.nq] + background 
     400    __call__ = Iq 
     401 
     402    def beta(self, call_details, values, cutoff, magnetic): 
     403        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
     404        self._call_kernel(call_details, values, cutoff, magnetic) 
     405        w_norm = self.result[2*self.q_input.nq + 1] 
     406        pd_norm = self.result[self.q_input.nq] 
     407        if w_norm == 0.: 
     408            w_norm = 1. 
     409        F2 = self.result[:self.q_input.nq]/w_norm 
     410        F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm 
     411        volume_avg = pd_norm/w_norm 
     412        return F1, F2, volume_avg 
     413 
     414    def _call_kernel(self, call_details, values, cutoff, magnetic): 
    390415        kernel = self.kernel[1 if magnetic else 0] 
    391416        args = [ 
    … …  
    394419            None, # pd_stop pd_stride[MAX_PD] 
    395420            call_details.buffer.ctypes.data, # problem 
    396             values.ctypes.data,  #pars 
    397             self.q_input.q.ctypes.data, #q 
     421            values.ctypes.data,  # pars 
     422            self.q_input.q.ctypes.data, # q 
    398423            self.result.ctypes.data,   # results 
    399424            self.real(cutoff), # cutoff 
    … …  
    407432            kernel(*args) # type: ignore 
    408433 
    409         #print("returned",self.q_input.q, self.result) 
    410         pd_norm = self.result[self.q_input.nq] 
    411         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    412         background = values[1] 
    413         #print("scale",scale,background) 
    414         return scale*self.result[:self.q_input.nq] + background 
    415  
    416434    def release(self): 
    417435        # type: () -> None 

  • sasmodels/modelinfo.py

    r1711569 rc036ddb  
    163163    parameter.length = length 
    164164    parameter.length_control = control 
    165  
    166165    return parameter 
    167166 
    … …  
    423422        self.kernel_parameters = parameters 
    424423        self._set_vector_lengths() 
    425  
    426424        self.npars = sum(p.length for p in self.kernel_parameters) 
    427425        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
    … …  
    430428        if self.nmagnetic: 
    431429            self.nvalues += 3 + 3*self.nmagnetic 
    432  
    433430        self.call_parameters = self._get_call_parameters() 
    434431        self.defaults = self._get_defaults() 
    … …  
    466463        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    467464                                if p.id.startswith('M0:')] 
    468  
    469465        self.pd_1d = set(p.name for p in self.call_parameters 
    470466                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    … …  
    770766        # Custom sum/multi models 
    771767        return kernel_module.model_info 
     768 
    772769    info = ModelInfo() 
    773770    #print("make parameter table", kernel_module.parameters) 
    … …  
    791788    info.category = getattr(kernel_module, 'category', None) 
    792789    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
     790    # TODO: find Fq by inspection 
     791    info.have_Fq = getattr(kernel_module, 'have_Fq', False) 
    793792    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    794793    info.source = getattr(kernel_module, 'source', []) 
    … …  
    822821    info.lineno = {} 
    823822    _find_source_lines(info, kernel_module) 
    824  
    825823    return info 
    826824 
    … …  
    915913    #: provided in the file. 
    916914    structure_factor = None # type: bool 
     915    #: True if the model defines an Fq function with signature 
     916    #: void Fq(double q, double *F1, double *F2, ...) 
     917    have_Fq = False 
    917918    #: List of C source files used to define the model.  The source files 
    918919    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 

  • sasmodels/models/core_shell_sphere.py

    rdc76240 rc036ddb  
    5858title = "Form factor for a monodisperse spherical particle with particle with a core-shell structure." 
    5959description = """ 
    60     F^2(q) = 3/V_s [V_c (sld_core-sld_shell) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
    61                    + V_s (sld_shell-sld_solvent) (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3] 
     60    F(q) = [V_c (sld_core-sld_shell) 3 (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
     61            + V_s (sld_shell-sld_solvent) 3 (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3] 
    6262 
    6363            V_s: Volume of the sphere shell 

  • sasmodels/models/ellipsoid.c

    r108e70e rc036ddb  
    55} 
    66 
     7/* Fq overrides Iq 
    78static  double 
    89Iq(double q, 
    … …  
    3738    return 1.0e-4 * s * s * form; 
    3839} 
     40*/ 
     41 
     42static void 
     43Fq(double q, 
     44    double *F1, 
     45    double *F2, 
     46    double sld, 
     47    double sld_solvent, 
     48    double radius_polar, 
     49    double radius_equatorial) 
     50{ 
     51    // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955) 
     52    //     i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT 
     53    //          = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT 
     54    //          = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT 
     55    // u-substitution of 
     56    //     u = sin, du = cos dT 
     57    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
     58    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
     59    // translate a point in [-1,1] to a point in [0, 1] 
     60    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     61    const double zm = 0.5; 
     62    const double zb = 0.5; 
     63    double total_F2 = 0.0; 
     64    double total_F1 = 0.0; 
     65    for (int i=0;i<GAUSS_N;i++) { 
     66        const double u = GAUSS_Z[i]*zm + zb; 
     67        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
     68        const double f = sas_3j1x_x(q*r); 
     69        total_F2 += GAUSS_W[i] * f * f; 
     70        total_F1 += GAUSS_W[i] * f; 
     71    } 
     72    // translate dx in [-1,1] to dx in [lower,upper] 
     73    const double form_squared_avg = total_F2*zm; 
     74    const double form_avg = total_F1*zm; 
     75    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
     76    *F1 = 1e-2 * s * form_avg; 
     77    *F2 = 1e-4 * s * s * form_squared_avg; 
     78} 
    3979 
    4080static double 

  • sasmodels/models/ellipsoid.py

    r2d81cfe rc036ddb  
    161161             ] 
    162162 
     163 
    163164source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
     165 
     166have_Fq = True 
    164167 
    165168def ER(radius_polar, radius_equatorial): 

  • sasmodels/models/lib/sphere_form.c

    r925ad6e r01c8d9e  
    1313    return 1.0e-4*square(contrast * fq); 
    1414} 
     15 

  • sasmodels/models/sphere.py

    ref07e95 rc036ddb  
    6969source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c"] 
    7070 
    71 # No volume normalization despite having a volume parameter 
    72 # This should perhaps be volume normalized? 
    73 form_volume = """ 
     71c_code = """ 
     72static double form_volume(double radius) 
     73{ 
    7474    return sphere_volume(radius); 
    75     """ 
     75} 
    7676 
    77 Iq = """ 
    78     return sphere_form(q, radius, sld, sld_solvent); 
    79     """ 
     77static void Fq(double q, double *F1,double *F2, double sld, double solvent_sld, double radius) 
     78{ 
     79    const double fq = sas_3j1x_x(q*radius); 
     80    const double contrast = (sld - solvent_sld); 
     81    const double form = 1e-2 * contrast * sphere_volume(radius) * fq; 
     82    *F1 = form; 
     83    *F2 = form*form; 
     84} 
     85""" 
     86 
     87# TODO: figure this out by inspection 
     88have_Fq = True 
    8089 
    8190def ER(radius): 

  • sasmodels/product.py

    r2d81cfe rc036ddb  
    1616import numpy as np  # type: ignore 
    1717 
    18 from .modelinfo import ParameterTable, ModelInfo 
     18from .modelinfo import ParameterTable, ModelInfo, Parameter 
    1919from .kernel import KernelModel, Kernel 
    2020from .details import make_details, dispersion_mesh 
    … …  
    3737ER_ID = "radius_effective" 
    3838VF_ID = "volfraction" 
     39BETA_DEFINITION = ("beta_mode", "", 0, [["P*S"],["P*(1+beta*(S-1))"]], "", 
     40                   "Structure factor dispersion calculation mode") 
    3941 
    4042# TODO: core_shell_sphere model has suppressed the volume ratio calculation 
    … …  
    7476    translate_name = dict((old.id, new.id) for old, new 
    7577                          in zip(s_pars.kernel_parameters[1:], s_list)) 
    76     combined_pars = p_pars.kernel_parameters + s_list 
     78    beta = [Parameter(*BETA_DEFINITION)] if p_info.have_Fq else [] 
     79    combined_pars = p_pars.kernel_parameters + s_list + beta 
    7780    parameters = ParameterTable(combined_pars) 
    7881    parameters.max_pd = p_pars.max_pd + s_pars.max_pd 
    … …  
    151154        #: Structure factor modelling interaction between particles. 
    152155        self.S = S 
     156 
    153157        #: Model precision. This is not really relevant, since it is the 
    154158        #: individual P and S models that control the effective dtype, 
    … …  
    168172        # in opencl; or both in opencl, but one in single precision and the 
    169173        # other in double precision). 
     174 
    170175        p_kernel = self.P.make_kernel(q_vectors) 
    171176        s_kernel = self.S.make_kernel(q_vectors) 
    … …  
    211216        p_offset = call_details.offset[:p_npars] 
    212217        p_details = make_details(p_info, p_length, p_offset, nweights) 
     218 
    213219        # Set p scale to the volume fraction in s, which is the first of the 
    214220        # 'S' parameters in the parameter list, or 2+np in 0-origin. 
    … …  
    254260        s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 
    255261 
     262        # beta mode is the first parameter after the structure factor pars 
     263        beta_index = 2+p_npars+s_npars 
     264        beta_mode = values[beta_index] 
     265 
    256266        # Call the kernels 
    257         p_result = self.p_kernel(p_details, p_values, cutoff, magnetic) 
    258         s_result = self.s_kernel(s_details, s_values, cutoff, False) 
    259  
    260         #print("p_npars",p_npars,s_npars,p_er,s_vr,values[2+p_npars+1:2+p_npars+s_npars]) 
     267        s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
     268        scale, background = values[0], values[1] 
     269        if beta_mode: 
     270            F1, F2, volume_avg = self.p_kernel.beta(p_details, p_values, cutoff, magnetic) 
     271            combined_scale = scale*volfrac/volume_avg 
     272            # Define lazy results based on intermediate values. 
     273            # The return value for the calculation should be an ordered 
     274            # dictionary containing any result the user might want to see 
     275            # at the end of the calculation, including scalars, strings, 
     276            # and plottable data.  Don't want to build this structure during 
     277            # fits, only when displaying the final result (or a one-off 
     278            # computation which asks for it). 
     279            # Do not use the current hack of storing the intermediate values 
     280            # in self.results since that leads to awkward threading issues. 
     281            # Instead return the function along with the bundle of inputs. 
     282            # P and Q may themselves have intermediate results they want to 
     283            # include, such as A and B if P = A + B.  Might use this mechanism 
     284            # to return the computed effective radius as well. 
     285            #def lazy_results(Q, S, F1, F2, scale): 
     286            #    """ 
     287            #    beta = F1**2 / F2  # what about divide by zero errors? 
     288            #    return { 
     289            #        'P' : Data1D(Q, scale*F2), 
     290            #        'beta': Data1D(Q, beta), 
     291            #        'S' : Data1D(Q, S), 
     292            #        'Seff': Data1D(Q, 1 + beta*(S-1)), 
     293            #        'I' : Data1D(Q, scale*(F2 + (F1**2)*(S-1)) + background), 
     294            #    } 
     295            #lazy_pars = s_result, F1, F2, combined_scale 
     296            self.results = [F2, s_result] 
     297            final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background 
     298        else: 
     299            p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic) 
     300            # remember the parts for plotting later 
     301            self.results = [p_result, s_result] 
     302            final_result = scale*(p_result*s_result) + background 
     303 
    261304        #call_details.show(values) 
    262305        #print("values", values) 
    … …  
    265308        #s_details.show(s_values) 
    266309        #print("=>", s_result) 
    267  
    268         # remember the parts for plotting later 
    269         self.results = [p_result, s_result] 
    270  
    271310        #import pylab as plt 
    272311        #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') 
    273312        #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') 
    274313        #plt.figure() 
    275  
    276         return values[0]*(p_result*s_result) + values[1] 
     314        return final_result 
    277315 
    278316    def release(self): 
Note: See TracChangeset for help on using the changeset viewer.