Changeset c036ddb in sasmodels


Ignore:
Timestamp:
Aug 7, 2018 8:45:45 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
7e923c2
Parents:
7b0abf8
Message:

refactor so Iq is not needed if Fq is defined

Location:
sasmodels
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/core.py

    r01c8d9e rc036ddb  
    140140    used with functions in generate to build the docs or extract model info. 
    141141    """ 
    142  
    143142    if "+" in model_string: 
    144143        parts = [load_model_info(part) 
    … …  
    290289    else: 
    291290        numpy_dtype = np.dtype(dtype) 
     291 
    292292    # Make sure that the type is supported by opencl, otherwise use dll 
    293293    if platform == "ocl": 

  • sasmodels/data.py

    r01c8d9e rc036ddb  
    522522                  else 'log') 
    523523        plt.xscale(xscale) 
    524          
    525524        plt.xlabel("$q$/A$^{-1}$") 
    526525        plt.yscale(yscale) 

  • sasmodels/generate.py

    r01c8d9e rc036ddb  
    671671 
    672672 
    673 # type in IQXY pattern could be single, float, double, long double, ... 
    674673_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]", 
    675674                           flags=re.MULTILINE) 
    … …  
    701700    return 'qa' 
    702701 
    703 # type in IQXY pattern could be single, float, double, long double, ... 
     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. 
    704705_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE) 
    705 def has_Fq(source): 
     706def contains_Fq(source): 
     707    # type: (List[str]) -> bool 
     708    """ 
     709    Return True if C source defines "void Fq(". 
     710    """ 
    706711    for code in source: 
    707712        m = _FQ_PATTERN.search(code) 
    … …  
    739744    # for computing volume even if we allow non-disperse volume parameters. 
    740745    partable = model_info.parameters 
     746 
    741747    # Load templates and user code 
    742748    kernel_header = load_template('kernel_header.c') 
    743749    kernel_code = load_template('kernel_iq.c') 
    744750    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
     751 
    745752    # Build initial sources 
    746753    source = [] 
    … …  
    755762    q, qx, qy, qab, qa, qb, qc \ 
    756763        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
     764 
    757765    # Generate form_volume function, etc. from body only 
    758766    if isinstance(model_info.form_volume, str): 
    … …  
    804812    source.append(call_volume) 
    805813    model_refs = _call_pars("_v.", partable.iq_parameters) 
    806     #create varaible BETA to turn on and off beta approximation 
    807     BETA = has_Fq(source) 
    808     if not BETA: 
     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: 
    809820        pars = ",".join(["_q"] + model_refs) 
    810821        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
    811     else: 
    812         pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs) 
    813         call_iq = "#define CALL_IQ(_q, _F1, _F2, _v) Fq(%s)" % pars 
     822        clear_iq = "#undef CALL_IQ" 
    814823    if xy_mode == 'qabc': 
    815824        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    … …  
    820829        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    821830        clear_iqxy = "#undef CALL_IQ_AC" 
    822     elif xy_mode == 'qa': 
     831    elif xy_mode == 'qa' and not model_info.have_Fq: 
    823832        pars = ",".join(["_qa"] + model_refs) 
    824833        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    825834        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" 
    826843    elif xy_mode == 'qxy': 
    827844        orientation_refs = _call_pars("_v.", partable.orientation_parameters) 
    … …  
    840857               if p.type == 'sld'] 
    841858    # Fill in definitions for numbers of parameters 
    842     source.append("#define BETA %d" %(1 if BETA else 0)) 
    843859    source.append("#define MAX_PD %s"%partable.max_pd) 
    844860    source.append("#define NUM_PARS %d"%partable.npars) 
    … …  
    848864    source.append("#define PROJECTION %d"%PROJECTION) 
    849865    # TODO: allow mixed python/opencl kernels? 
    850     ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    851     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) 
    852868 
    853869    result = { 
    … …  
    855871        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    856872    } 
    857     #print(result['dll']) 
    858873    return result 
    859874 
    860875 
    861 def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
     876def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name): 
    862877    # type: ([str,str], str, str, str) -> List[str] 
    863878    code = kernel[0] 
    … …  
    869884        '#line 1 "%s Iq"' % path, 
    870885        code, 
    871         "#undef CALL_IQ", 
     886        clear_iq, 
    872887        "#undef KERNEL_NAME", 
    873888        ] 
    … …  
    10751090        model_info = make_model_info(kernel_module) 
    10761091        source = make_source(model_info) 
    1077         #print(source['dll']) 
     1092        print(source['dll']) 
    10781093 
    10791094 

  • sasmodels/kernel_iq.c

    r01c8d9e 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 
    … …  
    8587  out_spin = clip(out_spin, 0.0, 1.0); 
    8688  // Previous version of this function took the square root of the weights, 
    87   // under the assumption that  
     89  // under the assumption that 
    8890  // 
    8991  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
    … …  
    272274 
    273275// ==================== KERNEL CODE ======================== 
     276#define COMPUTE_F1_F2 defined(CALL_FQ) 
    274277kernel 
    275278void KERNEL_NAME( 
    … …  
    338341  // The code differs slightly between opencl and dll since opencl is only 
    339342  // seeing one q value (stored in the variable "this_result") while the dll 
    340   // version must loop over all q 
    341   #if BETA 
    342   double *beta_result = &(result[nq+1]); // = result + nq + 1 
    343   double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
    344   #endif 
     343  // version must loop over all q. 
    345344  #ifdef USE_OPENCL 
    346     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    347     double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
    348     #if BETA 
    349       double this_beta_result = (pd_start == 0 ? 0.0 : result[nq + 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 
    350354  #else // !USE_OPENCL 
    351     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 
    352361    if (pd_start == 0) { 
    353362      #ifdef USE_OPENMP 
    354363      #pragma omp parallel for 
    355364      #endif 
    356       for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
    357       #if BETA 
    358       for (int q_index=0; q_index < nq; q_index++) beta_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; 
    359369      #endif 
    360370    } 
    361371    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    362372#endif // !USE_OPENCL 
    363  
    364  
    365  
    366373 
    367374 
    … …  
    454461 
    455462 
    456 #if defined(CALL_IQ) 
    457   // unoriented 1D 
     463#if defined(CALL_FQ) // COMPUTE_F1_F2 is true 
     464  // unoriented 1D returning <F> and <F^2> 
    458465  double qk; 
    459   #if BETA == 0 
    460     #define FETCH_Q() do { qk = q[q_index]; } while (0) 
    461     #define BUILD_ROTATION() do {} while(0) 
    462     #define APPLY_ROTATION() do {} while(0) 
    463     #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 
    464  
    465   // unoriented 1D Beta 
    466   #elif BETA == 1 
    467     double F1, F2; 
    468     #define FETCH_Q() do { qk = q[q_index]; } while (0) 
    469     #define BUILD_ROTATION() do {} while(0) 
    470     #define APPLY_ROTATION() do {} while(0) 
    471     #define CALL_KERNEL() CALL_IQ(qk,F1,F2,local_values.table) 
    472   #endif 
     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> 
     483  double qk; 
     484  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     485  #define BUILD_ROTATION() do {} while(0) 
     486  #define APPLY_ROTATION() do {} while(0) 
     487  #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 
    473488 
    474489#elif defined(CALL_IQ_A) 
    … …  
    504519  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
    505520  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     521 
    506522#elif defined(CALL_IQ_XY) 
    507523  // direct call to qx,qy calculator 
    … …  
    517533// logic in the IQ_AC and IQ_ABC branches.  This will become more important 
    518534// if we implement more projections, or more complicated projections. 
    519 #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 
    520537  #define APPLY_PROJECTION() const double weight=weight0 
    521538#elif defined(CALL_IQ_XY) // pass orientation to the model 
    … …  
    668685    if (weight > cutoff) { 
    669686      pd_norm += weight * CALL_VOLUME(local_values.table); 
     687      #if COMPUTE_F1_F2 
     688      weight_norm += weight; 
     689      #endif 
    670690      BUILD_ROTATION(); 
    671 #if BETA 
    672     if (weight > cutoff) { 
    673       weight_norm += weight;} 
    674 #endif 
     691 
    675692#ifndef USE_OPENCL 
    676693      // DLL needs to explicitly loop over the q values. 
    … …  
    718735          } 
    719736        #else  // !MAGNETIC 
    720           #if defined(CALL_IQ) && BETA == 1 
    721             CALL_KERNEL(); 
    722             const double scatteringF1 = F1; 
    723             const double scatteringF2 = F2; 
     737          #if COMPUTE_F1_F2 
     738            CALL_KERNEL(); // sets F1 and F2 by reference 
    724739          #else 
    725740            const double scattering = CALL_KERNEL(); 
    … …  
    729744 
    730745        #ifdef USE_OPENCL 
    731           #if defined(CALL_IQ)&& BETA == 1 
    732              this_result += weight * scatteringF2; 
    733              this_beta_result += weight * scatteringF1; 
    734             #else 
    735               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; 
    736751          #endif 
    737752        #else // !USE_OPENCL 
    738           #if defined(CALL_IQ)&& BETA == 1 
    739             result[q_index] += weight * scatteringF2; 
    740             beta_result[q_index] += weight * scatteringF1; 
    741             #endif 
    742             #else 
     753          #if COMPUTE_F1_F2 
     754            result[q_index] += weight * F2; 
     755            result[q_index+nq+1] += weight * F1; 
     756          #else 
    743757            result[q_index] += weight * scattering; 
    744758          #endif 
    … …  
    767781// Remember the current result and the updated norm. 
    768782#ifdef USE_OPENCL 
    769   result[q_index] = this_result; 
    770   if (q_index == 0) result[nq] = pd_norm; 
    771   #if BETA 
    772   beta_result[q_index] = this_beta_result; 
     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; 
    773793  #endif 
    774   if (q_index == 0) result[nq] = pd_norm; 
    775794 
    776795//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
    777796#else // !USE_OPENCL 
    778   result[nq] = pd_norm; 
    779   #if BETA 
    780   result[2*nq+1] = weight_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; 
    781802  #endif 
    782803//printf("res: %g/%g\n", result[0], pd_norm); 
    … …  
    784805 
    785806// ** clear the macros in preparation for the next kernel ** 
     807#undef COMPUTE_F1_F2 
    786808#undef PD_INIT 
    787809#undef PD_OPEN 

  • sasmodels/kernelcl.py

    r01c8d9e rc036ddb  
    420420        if self.program is None: 
    421421            compile_program = environment().compile_program 
    422             with open('model.c','w') as fid: 
    423                 print(self.source['opencl'], file=fid) 
    424422            timestamp = generate.ocl_timestamp(self.info) 
    425423            self.program = compile_program( 
    … …  
    540538        self.dtype = dtype 
    541539        self.dim = '2d' if q_input.is_2d else '1d' 
    542         # plus three for the normalization values 
    543         self.result = np.empty(2*q_input.nq+2,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) 
    544544 
    545545        # Inputs and outputs for each kernel call 
    … …  
    549549 
    550550        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    551                                   q_input.global_size[0] * dtype.itemsize) 
     551                                  q_input.global_size[0] * num_returns * dtype.itemsize) 
    552552        self.q_input = q_input # allocated by GpuInput above 
    553553 
    … …  
    557557                     else np.float16 if dtype == generate.F16 
    558558                     else np.float32)  # will never get here, so use np.float32 
    559     __call__= Iq 
    560559 
    561560    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): 
    562586        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    563587        context = self.queue.context 
    … …  
    575599        ] 
    576600        #print("Calling OpenCL") 
    577         call_details.show(values) 
     601        #call_details.show(values) 
    578602        #Call kernel and retrieve results 
    579603        wait_for = None 
    … …  
    601625                v.release() 
    602626 
    603         pd_norm = self.result[self.q_input.nq] 
    604         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    605         background = values[1] 
    606         #print("scale",scale,values[0],self.result[self.q_input.nq],background) 
    607         return scale*self.result[:self.q_input.nq] + background 
    608         # return self.result[:self.q_input.nq] 
    609      #NEEDS TO BE FINISHED FOR OPENCL 
    610      def beta(): 
    611          return 0 
    612  
    613627    def release(self): 
    614628        # type: () -> None 

  • sasmodels/kerneldll.py

    r01c8d9e rc036ddb  
    258258        # Note: if there is a syntax error then compile raises an error 
    259259        # and the source file will not be deleted. 
    260         #os.unlink(filename) 
    261         print("saving compiled file in %r"%filename) 
     260        os.unlink(filename) 
     261        #print("saving compiled file in %r"%filename) 
    262262    return dll 
    263263 
    … …  
    377377        self.dtype = q_input.dtype 
    378378        self.dim = '2d' if q_input.is_2d else '1d' 
    379         self.result = np.empty(2*q_input.nq+2, q_input.dtype) 
     379        # leave room for f1/f2 results in case we need to compute beta for 1d models 
     380        num_returns = 1 if self.dim == '2d' else 2  # 
     381        # plus 1 for the normalization value 
     382        self.result = np.empty((q_input.nq+1)*num_returns, self.dtype) 
    380383        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    381384                     else np.float64 if self.q_input.dtype == generate.F64 
    … …  
    417420            self.real(cutoff), # cutoff 
    418421        ] 
    419         #print(self.beta_result.ctypes.data) 
    420 #        print("Calling DLL line 397") 
    421 #        print("values", values) 
     422        #print("Calling DLL") 
    422423        #call_details.show(values) 
    423424        step = 100 

  • sasmodels/modelinfo.py

    r01c8d9e rc036ddb  
    417417    # scale and background are implicit parameters 
    418418    COMMON = [Parameter(*p) for p in COMMON_PARAMETERS] 
     419 
    419420    def __init__(self, parameters): 
    420421        # type: (List[Parameter]) -> None 
    … …  
    440441        self.form_volume_parameters = [p for p in self.kernel_parameters 
    441442                                       if p.type == 'volume'] 
     443 
    442444        # Theta offset 
    443445        offset = 0 
    … …  
    786788    info.category = getattr(kernel_module, 'category', None) 
    787789    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) 
    788792    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    789793    info.source = getattr(kernel_module, 'source', []) 
    … …  
    909913    #: provided in the file. 
    910914    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 
    911918    #: List of C source files used to define the model.  The source files 
    912919    #: 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

    r01c8d9e rc036ddb  
    55} 
    66 
     7/* Fq overrides Iq 
    78static  double 
    89Iq(double q, 
    … …  
    2021    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
    2122    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
    22    
     23 
    2324    // translate a point in [-1,1] to a point in [0, 1] 
    2425    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    … …  
    3637    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    3738    return 1.0e-4 * s * s * form; 
    38 }  
     39} 
     40*/ 
    3941 
    4042static void 
    … …  
    7274    const double form_avg = total_F1*zm; 
    7375    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
     76    *F1 = 1e-2 * s * form_avg; 
    7477    *F2 = 1e-4 * s * s * form_squared_avg; 
    75     *F1 = 1e-2 * s * form_avg; 
    7678} 
    77  
    78  
    79  
    80  
    81  
    8279 
    8380static 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/sphere.c

    r01c8d9e rc036ddb  
    44} 
    55 
     6// TODO: remove Iq when Iqxy can use Fq directly 
    67static double Iq(double q, double sld, double sld_solvent, double radius) 
    78{ 
    … …  
    1112static void Fq(double q, double *F1,double *F2, double sld, double solvent_sld, double radius) 
    1213{ 
    13     const double fq = sphere_volume(radius) * sas_3j1x_x(q*radius); 
     14    const double fq = sas_3j1x_x(q*radius); 
    1415    const double contrast = (sld - solvent_sld); 
    15     const double form = 1e-2*contrast * fq; 
    16     *F1= form; 
     16    const double form = 1e-2 * contrast * sphere_volume(radius) * fq; 
     17    *F1 = form; 
    1718    *F2 = form*form; 
    1819} 

  • sasmodels/models/sphere.py

    r01c8d9e rc036ddb  
    6767             ] 
    6868 
    69 source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c", "sphere.c"] 
     69source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c"] 
    7070 
     71c_code = """ 
     72static double form_volume(double radius) 
     73{ 
     74    return sphere_volume(radius); 
     75} 
     76 
     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 
    7189 
    7290def ER(radius): 

  • sasmodels/product.py

    r01c8d9e rc036ddb  
    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     beta_parameter = Parameter("beta_mode", "", 0, [["P*S"],["P*(1+beta*(S-1))"], "", "Structure factor dispersion calculation mode"]) 
    77     combined_pars = p_pars.kernel_parameters + s_list + [beta_parameter] 
     78    beta = [Parameter(*BETA_DEFINITION)] if p_info.have_Fq else [] 
     79    combined_pars = p_pars.kernel_parameters + s_list + beta 
    7880    parameters = ParameterTable(combined_pars) 
    7981    parameters.max_pd = p_pars.max_pd + s_pars.max_pd 
    … …  
    152154        #: Structure factor modelling interaction between particles. 
    153155        self.S = S 
    154          
     156 
    155157        #: Model precision. This is not really relevant, since it is the 
    156158        #: individual P and S models that control the effective dtype, 
    … …  
    170172        # in opencl; or both in opencl, but one in single precision and the 
    171173        # other in double precision). 
    172          
     174 
    173175        p_kernel = self.P.make_kernel(q_vectors) 
    174176        s_kernel = self.S.make_kernel(q_vectors) 
    … …  
    196198        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray 
    197199        p_info, s_info = self.info.composition[1] 
     200 
    198201        # if there are magnetic parameters, they will only be on the 
    199202        # form factor P, not the structure factor S. 
    … …  
    207210        nweights = call_details.num_weights 
    208211        weights = values[nvalues:nvalues + 2*nweights] 
     212 
    209213        # Construct the calling parameters for P. 
    210214        p_npars = p_info.parameters.npars 
    … …  
    212216        p_offset = call_details.offset[:p_npars] 
    213217        p_details = make_details(p_info, p_length, p_offset, nweights) 
     218 
    214219        # Set p scale to the volume fraction in s, which is the first of the 
    215220        # 'S' parameters in the parameter list, or 2+np in 0-origin. 
    … …  
    219224        p_values.append([0.]*spacer) 
    220225        p_values = np.hstack(p_values).astype(self.p_kernel.dtype) 
     226 
    221227        # Call ER and VR for P since these are needed for S. 
    222228        p_er, p_vr = calc_er_vr(p_info, p_details, p_values) 
    223229        s_vr = (volfrac/p_vr if p_vr != 0. else volfrac) 
    224230        #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr)) 
     231 
    225232        # Construct the calling parameters for S. 
    226233        # The  effective radius is not in the combined parameter list, so 
    … …  
    252259        s_values.append([0.]*spacer) 
    253260        s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 
     261 
    254262        # beta mode is the first parameter after the structure factor pars 
    255263        beta_index = 2+p_npars+s_npars 
    256264        beta_mode = values[beta_index] 
     265 
    257266        # Call the kernels 
    258         if beta_mode: # beta: 
     267        s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
     268        scale, background = values[0], values[1] 
     269        if beta_mode: 
    259270            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 
    260298        else: 
    261299            p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic) 
    262         s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    263         #print("p_npars",p_npars,s_npars,p_er,s_vr,values[2+p_npars+1:2+p_npars+s_npars]) 
     300            # remember the parts for plotting later 
     301            self.results = [p_result, s_result] 
     302            final_result = scale*(p_result*s_result) + background 
     303 
    264304        #call_details.show(values) 
    265305        #print("values", values) 
    … …  
    272312        #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') 
    273313        #plt.figure() 
    274         if beta_mode:#beta 
    275             beta_factor = F1**2/F2 
    276             Sq_eff = 1+beta_factor*(s_result - 1) 
    277             self.results = [F2, Sq_eff,F1,s_result] 
    278             final_result = volfrac*values[0]*(F2 + (F1**2)*(s_result - 1))/volume_avg+values[1] 
    279             #final_result =  volfrac * values[0] * F2 * Sq_eff / volume_avg + values[1] 
    280         else: 
    281             # remember the parts for plotting later 
    282             self.results = [p_result, s_result] 
    283             final_result = values[0]*(p_result*s_result) + values[1] 
    284314        return final_result 
    285315 
Note: See TracChangeset for help on using the changeset viewer.