Changeset 6e7ba14 in sasmodels


Ignore:
Timestamp:
Aug 20, 2018 6:52:21 AM (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:
aa44a6a
Parents:
bad3093
Message:

allow for different forms of effective_radius

Location:
sasmodels
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/generate.py

    rc036ddb r6e7ba14  
    805805        refs = _call_pars("_v.", partable.form_volume_parameters) 
    806806        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs)) 
     807        call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, _v) effective_radius(mode, %s)"%(",".join(refs)) 
    807808    else: 
    808809        # Model doesn't have volume.  We could make the kernel run a little 
    … …  
    810811        # the ifdef's reduce readability more than is worthwhile. 
    811812        call_volume = "#define CALL_VOLUME(v) 1.0" 
     813        call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, v) 0.0" 
    812814    source.append(call_volume) 
     815    source.append(call_effective_radius) 
    813816    model_refs = _call_pars("_v.", partable.iq_parameters) 
    814817 

  • sasmodels/kernel.py

    r2d81cfe r6e7ba14  
    4141    dtype = None  # type: np.dtype 
    4242 
    43     def __call__(self, call_details, values, cutoff, magnetic): 
     43    def Iq(self, call_details, values, cutoff, magnetic): 
    4444        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    45         raise NotImplementedError("need to implement __call__") 
     45        Pq, Reff = self.Pq_Reff(call_details, values, cutoff, magnetic, effective_radius_type=0) 
     46        scale = values[0] 
     47        background = values[1] 
     48        return scale*Pq + background 
     49    __call__ = Iq 
     50 
     51    def Pq_Reff(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     52        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     53        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     54        #print("returned",self.q_input.q, self.result) 
     55        nout = 2 if self.info.have_Fq else 1 
     56        total_weight = self.result[nout*self.q_input.nq + 0] 
     57        if total_weight == 0.: 
     58            total_weight = 1. 
     59        weighted_volume = self.result[nout*self.q_input.nq + 1] 
     60        weighted_radius = self.result[nout*self.q_input.nq + 2] 
     61        effective_radius = weighted_radius/total_weight 
     62        # compute I = scale*P + background 
     63        #           = scale*(sum(w*F^2)/sum w)/(sum (w*V)/sum w) + background 
     64        #           = scale/sum (w*V) * sum(w*F^2) + background 
     65        scale = values[0]/(weighted_volume if weighted_volume != 0.0 else 1.0) 
     66        background = values[1] 
     67        #print("scale",scale,background) 
     68        return self.result[0:nout*self.q_input.nq:nout], effective_radius 
     69 
     70    def beta(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     71        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     72        if self.dim == '2d': 
     73            raise NotImplementedError("beta not yet supported for 2D") 
     74        if not self.info.have_Fq: 
     75            raise NotImplementedError("beta not yet supported for "+self.info.id) 
     76        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     77        total_weight = self.result[2*self.q_input.nq + 0] 
     78        if total_weight == 0.: 
     79            total_weight = 1. 
     80        weighted_volume = self.result[2*self.q_input.nq + 1] 
     81        weighted_radius = self.result[2*self.q_input.nq + 2] 
     82        volume_average = weighted_volume/total_weight 
     83        effective_radius = weighted_radius/total_weight 
     84        F2 = self.result[0:2*self.q_input.nq:2]/total_weight 
     85        F1 = self.result[1:2*self.q_input.nq:2]/total_weight 
     86        return F1, F2, volume_average, effective_radius 
    4687 
    4788    def release(self): 
    4889        # type: () -> None 
    4990        pass 
     91 
     92    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     93        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     94        raise NotImpmentedError() 

  • sasmodels/kernel_iq.c

    rc036ddb r6e7ba14  
    2727//      parameters in the parameter table. 
    2828//  CALL_VOLUME(table) : call the form volume function 
     29//  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function 
    2930//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
    3031//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
    … …  
    284285    global const double *q, // nq q values, with padding to boundary 
    285286    global double *result,  // nq+1 return values, again with padding 
    286     const double cutoff     // cutoff in the dispersity weight product 
     287    const double cutoff,     // cutoff in the dispersity weight product 
     288    int effective_radius_type // which effective radius to compute 
    287289    ) 
    288290{ 
    … …  
    344346  #ifdef USE_OPENCL 
    345347    #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]); 
     348      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]); 
     349      double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     350      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]); 
     351      double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]); 
     352      double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]); 
    350353    #else 
    351       double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     354      double weight_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     355      double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]); 
     356      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]); 
    352357      double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
    353358    #endif 
    354359  #else // !USE_OPENCL 
    355360    #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]); 
     361      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]); 
     362      double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     363      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]); 
    358364    #else 
    359       double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     365      double weight_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     366      double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]); 
     367      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]); 
    360368    #endif 
    361369    if (pd_start == 0) { 
    … …  
    364372      #endif 
    365373      #if COMPUTE_F1_F2 
    366           for (int q_index=0; q_index < 2*nq+2; q_index++) result[q_index] = 0.0; 
     374          // 2*nq for F^2,F pairs 
     375          for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0; 
    367376      #else 
    368377          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
    369378      #endif 
    370379    } 
    371     //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
     380    //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_volume, result[0]); 
    372381#endif // !USE_OPENCL 
    373382 
    … …  
    684693    // Note: weight==0 must always be excluded 
    685694    if (weight > cutoff) { 
    686       pd_norm += weight * CALL_VOLUME(local_values.table); 
     695      weighted_volume += weight * CALL_VOLUME(local_values.table); 
    687696      #if COMPUTE_F1_F2 
    688697      weight_norm += weight; 
    689698      #endif 
     699      if (effective_radius_type != 0) { 
     700        weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); 
     701      } 
    690702      BUILD_ROTATION(); 
    691703 
    … …  
    745757        #ifdef USE_OPENCL 
    746758          #if COMPUTE_F1_F2 
     759            this_F2 += weight * F2; 
    747760            this_F1 += weight * F1; 
    748             this_F2 += weight * F2; 
    749761          #else 
    750762            this_result += weight * scattering; 
    … …  
    752764        #else // !USE_OPENCL 
    753765          #if COMPUTE_F1_F2 
    754             result[q_index] += weight * F2; 
    755             result[q_index+nq+1] += weight * F1; 
     766            result[2*q_index+0] += weight * F2; 
     767            result[2*q_index+1] += weight * F1; 
    756768          #else 
    757769            result[q_index] += weight * scattering; 
    … …  
    782794#ifdef USE_OPENCL 
    783795  #if COMPUTE_F1_F2 
    784     result[q_index] = this_F2; 
    785     result[q_index+nq+1] = this_F1; 
     796    result[2*q_index+0] = this_F2; 
     797    result[2*q_index+1] = this_F1; 
    786798    if (q_index == 0) { 
    787       result[nq] = pd_norm; 
    788       result[2*nq+1] = weight_norm; 
     799      result[2*nq+0] = weight_norm; 
     800      result[2*nq+1] = weighted_volume; 
     801      result[2*nq+2] = weighted_radius; 
    789802    } 
    790803  #else 
    791804    result[q_index] = this_result; 
    792     if (q_index == 0) result[nq] = pd_norm; 
     805    if (q_index == 0) { 
     806      result[nq+0] = weight_norm; 
     807      result[nq+1] = weighted_volume; 
     808      result[nq+2] = weighted_radius; 
     809    } 
    793810  #endif 
    794811 
    795 //if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
     812//if (q_index == 0) printf("res: %g/%g\n", result[0], weigthed_volume); 
    796813#else // !USE_OPENCL 
    797814  #if COMPUTE_F1_F2 
    798     result[nq] = pd_norm; 
    799     result[2*nq+1] = weight_norm; 
     815    result[2*nq] = weight_norm; 
     816    result[2*nq+1] = weighted_volume; 
     817    result[2*nq+2] = weighted_radius; 
    800818  #else 
    801     result[nq] = pd_norm; 
     819    result[nq] = weight_norm; 
     820    result[nq+1] = weighted_volume; 
     821    result[nq+2] = weighted_radius; 
    802822  #endif 
    803 //printf("res: %g/%g\n", result[0], pd_norm); 
     823//printf("res: %g/%g\n", result[0], weighted_volume); 
    804824#endif // !USE_OPENCL 
    805825 

  • sasmodels/kernelcl.py

    rc036ddb r6e7ba14  
    540540        # leave room for f1/f2 results in case we need to compute beta for 1d models 
    541541        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) 
     542        # plus 1 for the normalization value, plus another for R_eff 
     543        self.result = np.empty((q_input.nq+2)*num_returns, dtype) 
    544544 
    545545        # Inputs and outputs for each kernel call 
    … …  
    558558                     else np.float32)  # will never get here, so use np.float32 
    559559 
    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): 
     560    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    586561        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    587562        context = self.queue.context 
    … …  
    597572            details_b, values_b, self.q_input.q_b, self.result_b, 
    598573            self.real(cutoff), 
     574            np.uint32(effective_radius_type), 
    599575        ] 
    600576        #print("Calling OpenCL") 

  • sasmodels/kerneldll.py

    r2f8cbb9 r6e7ba14  
    315315 
    316316        # int, int, int, int*, double*, double*, double*, double*, double 
    317         argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type] 
     317        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type, ct.c_int32] 
    318318        names = [generate.kernel_name(self.info, variant) 
    319319                 for variant in ("Iq", "Iqxy", "Imagnetic")] 
    … …  
    382382        self.dim = '2d' if q_input.is_2d else '1d' 
    383383        # 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) 
     384        nout = 2 if self.info.have_Fq else 1 
     385        # plus 3 weight, volume, radius 
     386        self.result = np.empty(q_input.nq*nout + 3, self.dtype) 
    387387        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    388388                     else np.float64 if self.q_input.dtype == generate.F64 
    389389                     else np.float128) 
    390390 
    391     def Iq(self, call_details, values, cutoff, magnetic): 
    392         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    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         if self.dim == '2d': 
    405             raise NotImplementedError("beta not yet supported for 2D") 
    406         self._call_kernel(call_details, values, cutoff, magnetic) 
    407         w_norm = self.result[2*self.q_input.nq + 1] 
    408         pd_norm = self.result[self.q_input.nq] 
    409         if w_norm == 0.: 
    410             w_norm = 1. 
    411         F2 = self.result[:self.q_input.nq]/w_norm 
    412         F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm 
    413         volume_avg = pd_norm/w_norm 
    414         return F1, F2, volume_avg 
    415  
    416     def _call_kernel(self, call_details, values, cutoff, magnetic): 
     391    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     392        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    417393        kernel = self.kernel[1 if magnetic else 0] 
    418394        args = [ 
    … …  
    425401            self.result.ctypes.data,   # results 
    426402            self.real(cutoff), # cutoff 
     403            effective_radius_type, # cutoff 
    427404        ] 
    428405        #print("Calling DLL") 

  • sasmodels/kernelpy.py

    r108e70e r6e7ba14  
    1212 
    1313import numpy as np  # type: ignore 
     14from numpy import pi 
     15try: 
     16    from numpy import cbrt 
     17except ImportError: 
     18    def cbrt(x): return x ** (1.0/3.0) 
    1419 
    1520from .generate import F64 
    … …  
    158163 
    159164        # Generate a closure which calls the form_volume if it exists. 
    160         form_volume = model_info.form_volume 
    161         self._volume = ((lambda: form_volume(*volume_args)) if form_volume else 
    162                         (lambda: 1.0)) 
    163  
    164     def __call__(self, call_details, values, cutoff, magnetic): 
     165        self._volume_args = volume_args 
     166        volume = model_info.form_volume 
     167        radius = model_info.effective_radius 
     168        self._volume = ((lambda: volume(*volume_args)) if volume 
     169                        else (lambda: 1.0)) 
     170        self._radius = ((lambda mode: radius(mode, *volume_args)) if radius 
     171                        else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume 
     172                        else (lambda mode: 1.0)) 
     173 
     174 
     175 
     176    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    165177        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    166178        if magnetic: 
    … …  
    168180        #print("Calling python kernel") 
    169181        #call_details.show(values) 
    170         res = _loops(self._parameter_vector, self._form, self._volume, 
     182        radius = ((lambda: 0.0) if effective_radius_type == 0 
     183                  else (lambda: self._radius(effective_radius_type))) 
     184        res = _loops(self._parameter_vector, self._form, self._volume, radius, 
    171185                     self.q_input.nq, call_details, values, cutoff) 
    172186        return res 
    … …  
    183197           form,          # type: Callable[[], np.ndarray] 
    184198           form_volume,   # type: Callable[[], float] 
     199           form_radius,   # type: Callable[[], float] 
    185200           nq,            # type: int 
    186201           call_details,  # type: details.CallDetails 
    187202           values,        # type: np.ndarray 
    188            cutoff         # type: float 
     203           cutoff,        # type: float 
    189204          ): 
    190205    # type: (...) -> None 
    … …  
    209224    pd_weight = values[2+n_pars + call_details.num_weights:] 
    210225 
    211     pd_norm = 0.0 
     226    weight_norm = 0.0 
     227    weighted_volume = 0.0 
     228    weighted_radius = 0.0 
    212229    partial_weight = np.NaN 
    213230    weight = np.NaN 
    … …  
    246263            # update value and norm 
    247264            total += weight * Iq 
    248             pd_norm += weight * form_volume() 
    249  
    250     scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    251     background = values[1] 
    252     return scale*total + background 
     265            weight_norm += weight 
     266            weighted_volume += weight * form_volume() 
     267            weighted_radius += weight * form_radius() 
     268 
     269    result = np.hstack((total, weight_norm, weighted_volume, weighted_radius)) 
     270    return result 
    253271 
    254272 

  • sasmodels/modelinfo.py

    r7e923c2 r6e7ba14  
    789789    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
    790790    # TODO: find Fq by inspection 
     791    info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None) 
    791792    info.have_Fq = getattr(kernel_module, 'have_Fq', False) 
    792793    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    793794    info.source = getattr(kernel_module, 'source', []) 
    794795    info.c_code = getattr(kernel_module, 'c_code', None) 
     796    info.ER = None  # CRUFT 
     797    info.VR = None  # CRUFT 
    795798    # TODO: check the structure of the tests 
    796799    info.tests = getattr(kernel_module, 'tests', []) 
    797     info.ER = getattr(kernel_module, 'ER', None) # type: ignore 
    798     info.VR = getattr(kernel_module, 'VR', None) # type: ignore 
    799800    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 
    800801    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    … …  
    922923    #: use those functions.  Form factors are indicated by providing 
    923924    #: an :attr:`ER` function. 
     925    effective_radius_type = None   # type: List[str] 
     926    #: Returns the occupied volume and the total volume for each parameter set. 
     927    #: See :attr:`ER` for details on the parameters. 
    924928    source = None           # type: List[str] 
    925929    #: The set of tests that must pass.  The format of the tests is described 
    … …  
    942946    #: each parameter set.  Multiplicity parameters will be received as 
    943947    #: arrays, with one row per polydispersity condition. 
    944     ER = None               # type: Optional[Callable[[np.ndarray], np.ndarray]] 
    945     #: Returns the occupied volume and the total volume for each parameter set. 
    946     #: See :attr:`ER` for details on the parameters. 
    947     VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]] 
    948     #: Arbitrary C code containing supporting functions, etc., to be inserted 
    949     #: after everything in source.  This can include Iq and Iqxy functions with 
    950     #: the full function signature, including all parameters. 
    951948    c_code = None 
    952949    #: Returns the form volume for python-based models.  Form volume is needed 

  • sasmodels/models/ellipsoid.c

    r71b751d r6e7ba14  
    44    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    55} 
     6 
     7 
     8static double 
     9radius_from_volume(double radius_polar, double radius_equatorial) 
     10{ 
     11    return cbrt(radius_polar*radius_equatorial*radius_equatorial); 
     12} 
     13 
     14static double 
     15radius_from_curvature(double radius_polar, double radius_equatorial) 
     16{ 
     17    // Trivial cases 
     18    if (radius_polar == radius_equatorial) return radius_polar; 
     19    if (radius_polar * radius_equatorial == 0.)  return 0.; 
     20 
     21    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     22    const double ratio = (radius_polar < radius_equatorial  
     23                          ? radius_polar / radius_equatorial 
     24                          : radius_equatorial / radius_polar); 
     25    const double e1 = sqrt(1.0 - ratio*ratio); 
     26    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     27    const double bL = (1.0 + e1) / (1.0 - e1); 
     28    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     29    const double delta = 0.75 * b1 * b2; 
     30    const double ddd = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial * radius_equatorial; 
     31    return 0.5 * cbrt(ddd); 
     32} 
     33 
     34static double 
     35effective_radius(int mode, double radius_polar, double radius_equatorial) 
     36{ 
     37    if (mode == 1) { 
     38        return radius_from_curvature(radius_polar, radius_equatorial); 
     39    } else if (mode == 2) { 
     40        return radius_from_volume(radius_polar, radius_equatorial); 
     41    } else if (mode == 3) { 
     42        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial); 
     43    } else { 
     44        return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial); 
     45    } 
     46} 
     47 
    648 
    749static void 

  • sasmodels/models/ellipsoid.py

    r71b751d r6e7ba14  
    164164source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
    165165have_Fq = True 
    166  
    167 def ER(radius_polar, radius_equatorial): 
    168     # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
    169     ee = np.empty_like(radius_polar) 
    170     idx = radius_polar > radius_equatorial 
    171     ee[idx] = (radius_polar[idx] ** 2 - radius_equatorial[idx] ** 2) / radius_polar[idx] ** 2 
    172     idx = radius_polar < radius_equatorial 
    173     ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2 
    174     idx = radius_polar == radius_equatorial 
    175     ee[idx] = 2 * radius_polar[idx] 
    176     valid = (radius_polar * radius_equatorial != 0) 
    177     bd = 1.0 - ee[valid] 
    178     e1 = np.sqrt(ee[valid]) 
    179     b1 = 1.0 + np.arcsin(e1) / (e1 * np.sqrt(bd)) 
    180     bL = (1.0 + e1) / (1.0 - e1) 
    181     b2 = 1.0 + bd / 2 / e1 * np.log(bL) 
    182     delta = 0.75 * b1 * b2 
    183  
    184     ddd = np.zeros_like(radius_polar) 
    185     ddd[valid] = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial ** 2 
    186     return 0.5 * ddd ** (1.0 / 3.0) 
     166effective_radius_type = ["average curvature", "equivalent sphere", "min radius", "max radius"] 
    187167 
    188168def random(): 

  • sasmodels/product.py

    r6d90684 r6e7ba14  
    3535#] 
    3636 
    37 ER_ID = "radius_effective" 
    38 VF_ID = "volfraction" 
    39 BETA_DEFINITION = ("beta_mode", "", 0, [["P*S"],["P*(1+beta*(S-1))"]], "", 
    40                    "Structure factor dispersion calculation mode") 
     37RADIUS_ID = "radius_effective" 
     38VOLFRAC_ID = "volfraction" 
     39def make_extra_pars(p_info): 
     40    pars = [] 
     41    if p_info.have_Fq: 
     42        par = Parameter("structure_factor_mode", "", 0, [["P*S","P*(1+beta*(S-1))"]], "", 
     43                        "Structure factor calculation") 
     44        pars.append(par) 
     45    if p_info.effective_radius_type is not None: 
     46        par = Parameter("radius_effective_mode", "", 0, 
     47                        [["unconstrained"] + p_info.effective_radius_type], 
     48                        "", "Effective radius calculation") 
     49        pars.append(par) 
     50    return pars 
    4151 
    4252# TODO: core_shell_sphere model has suppressed the volume ratio calculation 
    … …  
    5262    # have any magnetic parameters 
    5363    if not len(s_info.parameters.kernel_parameters) >= 2: 
    54         raise TypeError("S needs {} and {} as its first parameters".format(ER_ID, VF_ID)) 
    55     if not s_info.parameters.kernel_parameters[0].id == ER_ID: 
    56         raise TypeError("S needs {} as first parameter".format(ER_ID)) 
    57     if not s_info.parameters.kernel_parameters[1].id == VF_ID: 
    58         raise TypeError("S needs {} as second parameter".format(VF_ID)) 
     64        raise TypeError("S needs {} and {} as its first parameters".format(RADIUS_ID, VOLFRAC_ID)) 
     65    if not s_info.parameters.kernel_parameters[0].id == RADIUS_ID: 
     66        raise TypeError("S needs {} as first parameter".format(RADIUS_ID)) 
     67    if not s_info.parameters.kernel_parameters[1].id == VOLFRAC_ID: 
     68        raise TypeError("S needs {} as second parameter".format(VOLFRAC_ID)) 
    5969    if not s_info.parameters.magnetism_index == []: 
    6070        raise TypeError("S should not have SLD parameters") 
    … …  
    6272    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
    6373 
    64     # Create list of parameters for the combined model.  Skip the first 
    65     # parameter of S, which we verified above is effective radius.  If there 
     74    # Create list of parameters for the combined model.  If there 
    6675    # are any names in P that overlap with those in S, modify the name in S 
    6776    # to distinguish it. 
    6877    p_set = set(p.id for p in p_pars.kernel_parameters) 
    6978    s_list = [(_tag_parameter(par) if par.id in p_set else par) 
    70               for par in s_pars.kernel_parameters[1:]] 
     79              for par in s_pars.kernel_parameters] 
    7180    # Check if still a collision after renaming.  This could happen if for 
    7281    # example S has volfrac and P has both volfrac and volfrac_S. 
    … …  
    7483        raise TypeError("name collision: P has P.name and P.name_S while S has S.name") 
    7584 
     85    # make sure effective radius is not a polydisperse parameter in product 
     86    s_list[0] = copy(s_list[0]) 
     87    s_list[0].polydisperse = False 
     88 
    7689    translate_name = dict((old.id, new.id) for old, new 
    77                           in zip(s_pars.kernel_parameters[1:], s_list)) 
    78     beta = [Parameter(*BETA_DEFINITION)] if p_info.have_Fq else [] 
    79     combined_pars = p_pars.kernel_parameters + s_list + beta 
     90                          in zip(s_pars.kernel_parameters, s_list)) 
     91    combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info) 
    8092    parameters = ParameterTable(combined_pars) 
    8193    parameters.max_pd = p_pars.max_pd + s_pars.max_pd 
    8294    def random(): 
    8395        combined_pars = p_info.random() 
    84         s_names = set(par.id for par in s_pars.kernel_parameters[1:]) 
     96        s_names = set(par.id for par in s_pars.kernel_parameters) 
    8597        combined_pars.update((translate_name[k], v) 
    8698                             for k, v in s_info.random().items() 
    … …  
    225237        p_values = np.hstack(p_values).astype(self.p_kernel.dtype) 
    226238 
    227         # Call ER and VR for P since these are needed for S. 
    228         p_er, p_vr = calc_er_vr(p_info, p_details, p_values) 
    229         s_vr = (volfrac/p_vr if p_vr != 0. else volfrac) 
    230         #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr)) 
    231  
    232239        # Construct the calling parameters for S. 
    233         # The  effective radius is not in the combined parameter list, so 
    234         # the number of 'S' parameters is one less than expected.  The 
    235         # computed effective radius needs to be added into the weights 
    236         # vector, especially since it is a polydisperse parameter in the 
    237         # stand-alone structure factor models.  We will added it at the 
    238         # end so the remaining offsets don't need to change. 
    239         s_npars = s_info.parameters.npars-1 
     240        s_npars = s_info.parameters.npars 
    240241        s_length = call_details.length[p_npars:p_npars+s_npars] 
    241242        s_offset = call_details.offset[p_npars:p_npars+s_npars] 
    … …  
    243244        s_offset = np.hstack((nweights, s_offset)) 
    244245        s_details = make_details(s_info, s_length, s_offset, nweights+1) 
    245         v, w = weights[:nweights], weights[nweights:] 
    246246        s_values = [ 
    247             # scale=1, background=0, radius_effective=p_er, volfraction=s_vr 
    248             [1., 0., p_er, s_vr], 
    249             # structure factor parameters start after scale, background and 
    250             # all the form factor parameters.  Skip the volfraction parameter 
    251             # as well, since it is computed elsewhere, and go to the end of the 
    252             # parameter list. 
    253             values[2+p_npars+1:2+p_npars+s_npars], 
    254             # no magnetism parameters to include for S 
    255             # add er into the (value, weights) pairs 
    256             v, [p_er], w, [1.0] 
     247            # scale=1, background=0, 
     248            [1., 0.], 
     249            values[2+p_npars:2+p_npars+s_npars], 
     250            weights, 
    257251        ] 
    258252        spacer = (32 - sum(len(v) for v in s_values)%32)%32 
    … …  
    261255 
    262256        # 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] 
     257        extra_offset = 2+p_npars+s_npars 
     258        if p_info.have_Fq: 
     259            beta_mode = values[extra_offset] 
     260            extra_offset += 1 
     261        else: 
     262            beta_mode = 0 
     263        if p_info.effective_radius_type is not None: 
     264            effective_radius_type = values[extra_offset] 
     265            extra_offset += 1 
     266        else: 
     267            effective_radius_type = 0 
    265268 
    266269        # Call the kernels 
    267         s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    268270        scale, background = values[0], values[1] 
    269271        if beta_mode: 
    270             F1, F2, volume_avg = self.p_kernel.beta(p_details, p_values, cutoff, magnetic) 
     272            F1, F2, volume_avg, effective_radius = self.p_kernel.beta( 
     273                p_details, p_values, cutoff, magnetic, effective_radius_type) 
     274            if effective_radius_type > 0: 
     275                s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 
     276            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    271277            combined_scale = scale*volfrac/volume_avg 
    272278            # Define lazy results based on intermediate values. 
    … …  
    297303            final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background 
    298304        else: 
    299             p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic) 
     305            p_result, effective_radius = self.p_kernel.Pq_Reff( 
     306                p_details, p_values, cutoff, magnetic, effective_radius_type) 
     307            if effective_radius_type > 0: 
     308                s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 
     309            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    300310            # remember the parts for plotting later 
    301311            self.results = [p_result, s_result] 

  • sasmodels/sasview_model.py

    rd533590 r6e7ba14  
    272272    attrs['category'] = model_info.category 
    273273    attrs['is_structure_factor'] = model_info.structure_factor 
    274     attrs['is_form_factor'] = model_info.ER is not None 
     274    attrs['is_form_factor'] = model_info.effective_radius_type is not None 
    275275    attrs['is_multiplicity_model'] = variants[0] > 1 
    276276    attrs['multiplicity_info'] = variants 
Note: See TracChangeset for help on using the changeset viewer.