Changeset aa44a6a in sasmodels


Ignore:
Timestamp:
Aug 21, 2018 11:38:57 AM (6 years ago)
Author:
Torin Cooper-Bennun <torin.cooper-bennun@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
3d6526d
Parents:
6e7ba14 (diff), d32de68 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

merge from beta_approx_lazy_results; include ER as intermediate result; fixed-choices params use units, not limits, to store choices

Location:
sasmodels
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/product.py

    r6e7ba14 raa44a6a  
    1313from __future__ import print_function, division 
    1414 
     15from collections import OrderedDict 
     16 
    1517from copy import copy 
    1618import numpy as np  # type: ignore 
    … …  
    2224# pylint: disable=unused-import 
    2325try: 
    24     from typing import Tuple 
     26    from typing import Tuple, Callable, Union 
    2527except ImportError: 
    2628    pass 
    … …  
    4042    pars = [] 
    4143    if p_info.have_Fq: 
    42         par = Parameter("structure_factor_mode", "", 0, [["P*S","P*(1+beta*(S-1))"]], "", 
     44        par = Parameter("structure_factor_mode", ["P*S","P*(1+beta*(S-1))"], 0, None, "", 
    4345                        "Structure factor calculation") 
    4446        pars.append(par) 
    4547    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") 
     48        par = Parameter("radius_effective_mode", ["unconstrained"] + p_info.effective_radius_type, 0, None, "", 
     49                        "Effective radius calculation") 
    4950        pars.append(par) 
    5051    return pars 
    … …  
    157158    return par 
    158159 
     160def _intermediates(P, S, effective_radius): 
     161    # type: (np.ndarray, np.ndarray, float) -> OrderedDict[str, np.ndarray] 
     162    """ 
     163    Returns intermediate results for standard product (P(Q)*S(Q)) 
     164    """ 
     165    return OrderedDict(( 
     166        ("P(Q)", P), 
     167        ("S(Q)", S), 
     168        ("effective_radius", effective_radius), 
     169    )) 
     170 
     171def _intermediates_beta(F1, F2, S, scale, bg, effective_radius): 
     172    # type: (np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, float) -> OrderedDict[str, Union[np.ndarray, float]] 
     173    """ 
     174    Returns intermediate results for beta approximation-enabled product. The result may be an array or a float. 
     175    """ 
     176    # TODO: 1. include calculated Q vector 
     177    # TODO: 2. consider implications if there are intermediate results in P(Q) 
     178    return OrderedDict(( 
     179        ("P(Q)", scale*F2), 
     180        ("S(Q)", S), 
     181        ("beta(Q)", F1**2 / F2), 
     182        ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), 
     183        ("effective_radius", effective_radius), 
     184        # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), 
     185    )) 
     186 
    159187class ProductModel(KernelModel): 
    160188    def __init__(self, model_info, P, S): 
    … …  
    276304            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    277305            combined_scale = scale*volfrac/volume_avg 
    278             # Define lazy results based on intermediate values. 
    279             # The return value for the calculation should be an ordered 
    280             # dictionary containing any result the user might want to see 
    281             # at the end of the calculation, including scalars, strings, 
    282             # and plottable data.  Don't want to build this structure during 
    283             # fits, only when displaying the final result (or a one-off 
    284             # computation which asks for it). 
    285             # Do not use the current hack of storing the intermediate values 
    286             # in self.results since that leads to awkward threading issues. 
    287             # Instead return the function along with the bundle of inputs. 
    288             # P and Q may themselves have intermediate results they want to 
    289             # include, such as A and B if P = A + B.  Might use this mechanism 
    290             # to return the computed effective radius as well. 
    291             #def lazy_results(Q, S, F1, F2, scale): 
    292             #    """ 
    293             #    beta = F1**2 / F2  # what about divide by zero errors? 
    294             #    return { 
    295             #        'P' : Data1D(Q, scale*F2), 
    296             #        'beta': Data1D(Q, beta), 
    297             #        'S' : Data1D(Q, S), 
    298             #        'Seff': Data1D(Q, 1 + beta*(S-1)), 
    299             #        'I' : Data1D(Q, scale*(F2 + (F1**2)*(S-1)) + background), 
    300             #    } 
    301             #lazy_pars = s_result, F1, F2, combined_scale 
    302             self.results = [F2*volfrac/volume_avg, s_result] 
     306 
     307            self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background, effective_radius) 
    303308            final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background 
     309 
    304310        else: 
    305311            p_result, effective_radius = self.p_kernel.Pq_Reff( 
    … …  
    309315            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    310316            # remember the parts for plotting later 
    311             self.results = [p_result, s_result] 
     317            self.results = lambda: _intermediates(p_result, s_result, effective_radius) 
    312318            final_result = scale*(p_result*s_result) + background 
    313319 

  • sasmodels/sasview_model.py

    r6e7ba14 raa44a6a  
    381381                continue 
    382382            self.params[p.name] = p.default 
    383             self.details[p.id] = [p.units, p.limits[0], p.limits[1]] 
     383            if p.limits and type(p.limits) is list and len(p.limits) > 1: 
     384                self.details[p.id] = [p.units if not p.choices else p.choices, p.limits[0], p.limits[1]] 
     385            else: 
     386                self.details[p.id] = [p.units if not p.choices else p.choices, None, None] 
    384387            if p.polydisperse: 
    385388                self.details[p.id+".width"] = [ 
    … …  
    616619        # so that it returns a results object containing all the bits: 
    617620        #     the A, B, C, ... of the composition model (and any subcomponents?) 
    618         #     the P and S of the product model, 
     621        #     the P and S of the product model 
    619622        #     the combined model before resolution smearing, 
    620623        #     the sasmodel before sesans conversion, 
    … …  
    638641            with calculation_lock: 
    639642                self._calculate_Iq(qx) 
    640                 return self._intermediate_results 
     643                # for compatibility with sasview 4.3 
     644                results = self._intermediate_results() 
     645                return results["P(Q)"], results["S(Q)"] 
    641646        else: 
    642647            return None 

  • 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(): 
Note: See TracChangeset for help on using the changeset viewer.