Changeset e44432d in sasmodels


Ignore:
Timestamp:
Oct 19, 2018 5:46:26 PM (3 weeks ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
beta_approx, py3, ticket-1015-gpu-mem-error, ticket-1157, ticket-608-user-defined-weights, ticket_1156
Children:
be43e39
Parents:
2586ab72
Message:

support hollow models in structure factor calculations

Location:
sasmodels
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/generate.py

    r2773c66 re44432d  
    715715    return False 
    716716 
     717_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE) 
     718def contains_shell_volume(source): 
     719    # type: (List[str]) -> bool 
     720    """ 
     721    Return True if C source defines "void Fq(". 
     722    """ 
     723    for code in source: 
     724        m = _SHELL_VOLUME_PATTERN.search(code) 
     725        if m is not None: 
     726            return True 
     727    return False 
     728 
    717729def _add_source(source, code, path, lineno=1): 
    718730    """ 
     
    767779        pars = partable.form_volume_parameters 
    768780        source.append(_gen_fn(model_info, 'form_volume', pars)) 
     781    if isinstance(model_info.shell_volume, str): 
     782        pars = partable.form_volume_parameters 
     783        source.append(_gen_fn(model_info, 'shell_volume', pars)) 
    769784    if isinstance(model_info.Iq, str): 
    770785        pars = [q] + partable.iq_parameters 
     
    779794        pars = [qa, qb, qc] + partable.iq_parameters 
    780795        source.append(_gen_fn(model_info, 'Iqabc', pars)) 
     796 
     797    # Check for shell_volume in source 
     798    is_hollow = contains_shell_volume(source) 
    781799 
    782800    # What kind of 2D model do we need?  Is it consistent with the parameters? 
     
    802820                              for p in partable.kernel_parameters)) 
    803821    # Define the function calls 
    804     call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, v) 0.0" 
     822    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0" 
    805823    if partable.form_volume_parameters: 
    806824        refs = _call_pars("_v.", partable.form_volume_parameters) 
    807         call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs)) 
     825        if is_hollow: 
     826            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2) 
     827        else: 
     828            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs)) 
    808829        if model_info.effective_radius_type: 
    809             call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, _v) effective_radius(mode, %s)"%(",".join(refs)) 
     830            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs)) 
    810831    else: 
    811832        # Model doesn't have volume.  We could make the kernel run a little 
    812833        # faster by not using/transferring the volume normalizations, but 
    813834        # the ifdef's reduce readability more than is worthwhile. 
    814         call_volume = "#define CALL_VOLUME(v) 1.0" 
     835        call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)" 
    815836    source.append(call_volume) 
    816837    source.append(call_effective_radius) 
  • sasmodels/kernel.py

    ree60aa7 re44432d  
    4343    def Iq(self, call_details, values, cutoff, magnetic): 
    4444        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    45         Pq, Reff = self.Pq_Reff(call_details, values, cutoff, magnetic, effective_radius_type=0) 
    46         return Pq 
     45        r""" 
     46        Returns I(q) from the polydisperse average scattering. 
     47 
     48        .. math:: 
     49 
     50            I(q) = \text{scale} \cdot P(q) + \text{background} 
     51 
     52        With the correct choice of model and contrast, setting *scale* to 
     53        the volume fraction $V_f$ of particles should match the measured 
     54        absolute scattering.  Some models (e.g., vesicle) have volume fraction 
     55        built into the model, and do not need an additional scale. 
     56        """ 
     57        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic, 
     58                              effective_radius_type=0) 
     59        combined_scale = values[0]/shell_volume 
     60        background = values[1] 
     61        return combined_scale*F2 + background 
    4762    __call__ = Iq 
    4863 
    49     def Pq_Reff(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     64    def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0): 
    5065        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     66        r""" 
     67        Returns <F(q)>, <F(q)^2>, effective radius, shell volume and 
     68        form:shell volume ratio.  The <F(q)> term may be None if the 
     69        form factor does not support direct computation of $F(q)$ 
     70 
     71        $P(q) = <F^2(q)>/<V>$ is used for structure factor calculations, 
     72 
     73        .. math:: 
     74 
     75            I(q) = \text{scale} \cdot P(q) \cdot S(q) + \text{background} 
     76 
     77        For the beta approximation, this becomes 
     78 
     79        .. math:: 
     80 
     81            I(q) = \text{scale} * P (1 + <F>^2/<F^2> (S - 1)) + \text{background} 
     82                 = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background} 
     83 
     84        $<F(q)>$ and $<F^2(q)>$ are averaged by polydispersity in shape 
     85        and orientation, with each configuration $x_k$ having form factor 
     86        $F(q, x_k)$, weight $w_k$ and volume $V_k$.  The result is: 
     87 
     88        .. math:: 
     89 
     90            P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} 
     91 
     92        The form factor itself is scaled by volume and contrast to compute the 
     93        total scattering.  This is then squared, and the volume weighted 
     94        F^2 is then normalized by volume F.  For a given density, the number 
     95        of scattering centers is assumed to scale linearly with volume.  Later 
     96        scaling the resulting $P(q)$ by the volume fraction of particles 
     97        gives the total scattering on an absolute scale. Most models 
     98        incorporate the volume fraction into the overall scale parameter.  An 
     99        exception is vesicle, which includes the volume fraction parameter in 
     100        the model itself, scaling $F$ by $\surd V_f$ so that the math for 
     101        the beta approximation works out. 
     102 
     103        By scaling $P(q)$ by total weight $\sum w_k$, there is no need to make 
     104        sure that the polydisperisity distributions normalize to one.  In 
     105        particular, any distibution values $x_k$ outside the valid domain 
     106        of $F$ will not be included, and the distribution will be implicitly 
     107        truncated.  This is controlled by the parameter limits defined in the 
     108        model (which truncate the distribution before calling the kernel) as 
     109        well as any region excluded using the *INVALID* macro defined within 
     110        the model itself. 
     111 
     112        The volume used in the polydispersity calculation is the form volume 
     113        for solid objects or the shell volume for hollow objects.  Shell 
     114        volume should be used within $F$ so that the normalizing scale 
     115        represents the volume fraction of the shell rather than the entire 
     116        form.  This corresponds to the volume fraction of shell-forming 
     117        material added to the solvent. 
     118 
     119        The calculation of $S$ requires the effective radius and the 
     120        volume fraction of the particles.  The model can have several 
     121        different ways to compute effective radius, with the 
     122        *effective_radius_type* parameter used to select amongst them.  The 
     123        volume fraction of particles should be determined from the total 
     124        volume fraction of the form, not just the shell volume fraction. 
     125        This makes a difference for hollow shapes, which need to scale 
     126        the volume fraction by the returned volume ratio when computing $S$. 
     127        For solid objects, the shell volume is set to the form volume so 
     128        this scale factor evaluates to one and so can be used for both 
     129        hollow and solid shapes. 
     130        """ 
    51131        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
    52132        #print("returned",self.q_input.q, self.result) 
     
    55135        if total_weight == 0.: 
    56136            total_weight = 1. 
    57         weighted_volume = self.result[nout*self.q_input.nq + 1] 
    58         weighted_radius = self.result[nout*self.q_input.nq + 2] 
    59         effective_radius = weighted_radius/total_weight 
    60         # compute I = scale*P + background 
    61         #           = scale*(sum(w*F^2)/sum w)/(sum (w*V)/sum w) + background 
    62         #           = scale/sum (w*V) * sum(w*F^2) + background 
    63         F2 = self.result[0:nout*self.q_input.nq:nout] 
    64         scale = values[0]/(weighted_volume if weighted_volume != 0.0 else 1.0) 
    65         background = values[1] 
    66         Pq = scale*F2 + background 
    67         #print("scale",scale,background) 
    68         return Pq, 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 
     137        # Note: shell_volume == form_volume for solid objects 
     138        form_volume = self.result[nout*self.q_input.nq + 1]/total_weight 
     139        shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight 
     140        effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight 
     141        if shell_volume == 0.: 
     142            shell_volume = 1. 
     143        F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None 
     144        F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight 
     145        return F1, F2, effective_radius, shell_volume, form_volume/shell_volume 
    87146 
    88147    def release(self): 
     
    92151    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    93152        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     153        """ 
     154        Call the kernel.  Subclasses defining kernels for particular execution 
     155        engines need to provide an implementation for this. 
     156        """ 
    94157        raise NotImplementedError() 
  • sasmodels/kernel_iq.c

    r2773c66 re44432d  
    2626//  MAGNETIC_PARS : a comma-separated list of indices to the sld 
    2727//      parameters in the parameter table. 
    28 //  CALL_VOLUME(table) : call the form volume function 
     28//  CALL_VOLUME(form, shell, table) : assign form and shell values 
    2929//  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function 
    3030//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
     
    275275 
    276276// ==================== KERNEL CODE ======================== 
    277 #define COMPUTE_F1_F2 defined(CALL_FQ) 
    278277kernel 
    279278void KERNEL_NAME( 
     
    345344  // version must loop over all q. 
    346345  #ifdef USE_OPENCL 
    347     #if COMPUTE_F1_F2 
     346    #if defined(CALL_FQ) 
    348347      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]); 
     348      double weighted_form = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     349      double weighted_shell = (pd_start == 0 ? 0.0 : result[2*nq+2]); 
     350      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+3]); 
    351351      double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]); 
    352352      double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]); 
    353353    #else 
    354354      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]); 
     355      double weighted_form = (pd_start == 0 ? 0.0 : result[nq+1]); 
     356      double weighted_shell = (pd_start == 0 ? 0.0 : result[nq+2]); 
     357      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+3]); 
    357358      double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
    358359    #endif 
    359360  #else // !USE_OPENCL 
    360     #if COMPUTE_F1_F2 
     361    #if defined(CALL_FQ) 
    361362      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]); 
     363      double weighted_form = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     364      double weighted_shell = (pd_start == 0 ? 0.0 : result[2*nq+2]); 
     365      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+3]); 
    364366    #else 
    365367      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]); 
     368      double weighted_form = (pd_start == 0 ? 0.0 : result[nq+1]); 
     369      double weighted_shell = (pd_start == 0 ? 0.0 : result[nq+2]); 
     370      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+3]); 
    368371    #endif 
    369372    if (pd_start == 0) { 
     
    371374      #pragma omp parallel for 
    372375      #endif 
    373       #if COMPUTE_F1_F2 
     376      #if defined(CALL_FQ) 
    374377          // 2*nq for F^2,F pairs 
    375378          for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0; 
     
    378381      #endif 
    379382    } 
    380     //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_volume, result[0]); 
     383    //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_shell, result[0]); 
    381384#endif // !USE_OPENCL 
    382385 
     
    693696    // Note: weight==0 must always be excluded 
    694697    if (weight > cutoff) { 
    695       weighted_volume += weight * CALL_VOLUME(local_values.table); 
    696       #if COMPUTE_F1_F2 
     698      double form, shell; 
     699      CALL_VOLUME(form, shell, local_values.table); 
    697700      weight_norm += weight; 
    698       #endif 
     701      weighted_form += weight * form; 
     702      weighted_shell += weight * shell; 
    699703      if (effective_radius_type != 0) { 
    700704        weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); 
     
    747751          } 
    748752        #else  // !MAGNETIC 
    749           #if COMPUTE_F1_F2 
     753          #if defined(CALL_FQ) 
    750754            CALL_KERNEL(); // sets F1 and F2 by reference 
    751755          #else 
     
    756760 
    757761        #ifdef USE_OPENCL 
    758           #if COMPUTE_F1_F2 
     762          #if defined(CALL_FQ) 
    759763            this_F2 += weight * F2; 
    760764            this_F1 += weight * F1; 
     
    763767          #endif 
    764768        #else // !USE_OPENCL 
    765           #if COMPUTE_F1_F2 
     769          #if defined(CALL_FQ) 
    766770            result[2*q_index+0] += weight * F2; 
    767771            result[2*q_index+1] += weight * F1; 
     
    793797// Remember the current result and the updated norm. 
    794798#ifdef USE_OPENCL 
    795   #if COMPUTE_F1_F2 
     799  #if defined(CALL_FQ) 
    796800    result[2*q_index+0] = this_F2; 
    797801    result[2*q_index+1] = this_F1; 
    798802    if (q_index == 0) { 
    799803      result[2*nq+0] = weight_norm; 
    800       result[2*nq+1] = weighted_volume; 
    801       result[2*nq+2] = weighted_radius; 
     804      result[2*nq+1] = weighted_form; 
     805      result[2*nq+3] = weighted_shell; 
     806      result[2*nq+3] = weighted_radius; 
    802807    } 
    803808  #else 
     
    805810    if (q_index == 0) { 
    806811      result[nq+0] = weight_norm; 
    807       result[nq+1] = weighted_volume; 
    808       result[nq+2] = weighted_radius; 
     812      result[nq+1] = weighted_form; 
     813      result[nq+2] = weighted_shell; 
     814      result[nq+3] = weighted_radius; 
    809815    } 
    810816  #endif 
    811817 
    812 //if (q_index == 0) printf("res: %g/%g\n", result[0], weigthed_volume); 
     818//if (q_index == 0) printf("res: %g/%g\n", result[0], weighted_shell); 
    813819#else // !USE_OPENCL 
    814   #if COMPUTE_F1_F2 
     820  #if defined(CALL_FQ) 
    815821    result[2*nq] = weight_norm; 
    816     result[2*nq+1] = weighted_volume; 
    817     result[2*nq+2] = weighted_radius; 
     822    result[2*nq+1] = weighted_form; 
     823    result[2*nq+2] = weighted_shell; 
     824    result[2*nq+3] = weighted_radius; 
    818825  #else 
    819826    result[nq] = weight_norm; 
    820     result[nq+1] = weighted_volume; 
    821     result[nq+2] = weighted_radius; 
     827    result[nq+1] = weighted_form; 
     828    result[nq+2] = weighted_shell; 
     829    result[nq+3] = weighted_radius; 
    822830  #endif 
    823 //printf("res: %g/%g\n", result[0], weighted_volume); 
     831//printf("res: %g/%g\n", result[0], weighted_shell); 
    824832#endif // !USE_OPENCL 
    825833 
    826834// ** clear the macros in preparation for the next kernel ** 
    827 #undef COMPUTE_F1_F2 
    828835#undef PD_INIT 
    829836#undef PD_OPEN 
  • sasmodels/kernelcl.py

    r5399809 re44432d  
    539539        # leave room for f1/f2 results in case we need to compute beta for 1d models 
    540540        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
    541         # plus 3 weight, volume, radius 
    542         self.result = np.empty(q_input.nq*nout + 3, self.dtype) 
     541        # +4 for total weight, shell volume, effective radius, form volume 
     542        self.result = np.empty(q_input.nq*nout + 4, self.dtype) 
    543543 
    544544        # Inputs and outputs for each kernel call 
  • sasmodels/kerneldll.py

    r6e7ba14 re44432d  
    9999    pass 
    100100# pylint: enable=unused-import 
     101 
     102# Compiler output is a byte stream that needs to be decode in python 3 
     103decode = (lambda s: s) if sys.version_info[0] < 3 else (lambda s: s.decode('utf8')) 
    101104 
    102105if "SAS_DLL_PATH" in os.environ: 
     
    184187        subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT) 
    185188    except subprocess.CalledProcessError as exc: 
    186         raise RuntimeError("compile failed.\n%s\n%s"%(command_str, exc.output)) 
     189        output = decode(exc.output) 
     190        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, output)) 
    187191    if not os.path.exists(output): 
    188192        raise RuntimeError("compile failed.  File is in %r"%source) 
     
    383387        # leave room for f1/f2 results in case we need to compute beta for 1d models 
    384388        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) 
     389        # +4 for total weight, shell volume, effective radius, form volume 
     390        self.result = np.empty(q_input.nq*nout + 4, self.dtype) 
    387391        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    388392                     else np.float64 if self.q_input.dtype == generate.F64 
  • sasmodels/kernelpy.py

    r9ae0d2e re44432d  
    164164        self._volume_args = volume_args 
    165165        volume = model_info.form_volume 
     166        shell = model_info.shell_volume 
    166167        radius = model_info.effective_radius 
    167         self._volume = ((lambda: volume(*volume_args)) if volume 
    168                         else (lambda: 1.0)) 
     168        self._volume = ((lambda: (shell(*volume_args), volume(*volume_args))) if shell and volume 
     169                        else (lambda: [volume(*volume_args)]*2) if volume 
     170                        else (lambda: (1.0, 1.0))) 
    169171        self._radius = ((lambda mode: radius(mode, *volume_args)) if radius 
    170172                        else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume 
     
    225227        total = form() 
    226228        weight_norm = 1.0 
    227         weighted_volume = form_volume() 
     229        weighted_shell, weighted_form = form_volume() 
    228230        weighted_radius = form_radius() 
    229231 
     
    233235 
    234236        weight_norm = 0.0 
    235         weighted_volume = 0.0 
     237        weighted_form = 0.0 
     238        weighted_shell = 0.0 
    236239        weighted_radius = 0.0 
    237240        partial_weight = np.NaN 
     
    272275                total += weight * Iq 
    273276                weight_norm += weight 
    274                 weighted_volume += weight * form_volume() 
     277                shell, form = form_volume() 
     278                weighted_form += weight * form 
     279                weighted_shell += weight * shell 
    275280                weighted_radius += weight * form_radius() 
    276281 
    277     result = np.hstack((total, weight_norm, weighted_volume, weighted_radius)) 
     282    result = np.hstack((total, weight_norm, weighted_form, weighted_shell, weighted_radius)) 
    278283    return result 
    279284 
  • sasmodels/modelinfo.py

    r9ae0d2e re44432d  
    718718 
    719719#: Set of variables defined in the model that might contain C code 
    720 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
     720C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code'] 
    721721 
    722722def _find_source_lines(model_info, kernel_module): 
     
    803803    info.tests = getattr(kernel_module, 'tests', []) 
    804804    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 
     805    info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore 
    805806    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    806807    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 
     
    921922    #: void Fq(double q, double *F1, double *F2, ...) 
    922923    have_Fq = False 
     924    #: List of options for computing the effective radius of the shape, 
     925    #: or None if the model is not usable as a form factor model. 
     926    effective_radius_type = None   # type: List[str] 
    923927    #: List of C source files used to define the model.  The source files 
    924928    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 
    925929    #: model defines orientation parameters. Files containing the most basic 
    926930    #: functions must appear first in the list, followed by the files that 
    927     #: use those functions.  Form factors are indicated by providing 
    928     #: an :attr:`ER` function. 
    929     effective_radius_type = None   # type: List[str] 
    930     #: Returns the occupied volume and the total volume for each parameter set. 
    931     #: See :attr:`ER` for details on the parameters. 
     931    #: use those functions. 
    932932    source = None           # type: List[str] 
    933933    #: The set of tests that must pass.  The format of the tests is described 
     
    958958    #: C code, either defined as a string, or in the sources. 
    959959    form_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
     960    #: Returns the shell volume for python-based models.  Form volume and 
     961    #: shell volume are needed for volume normalization in the polydispersity 
     962    #: integral and structure interactions for hollow shapes.  If no 
     963    #: parameters are *volume* parameters, then shell volume is not needed. 
     964    #: For C-based models, (with :attr:`sources` defined, or with :attr:`Iq` 
     965    #: defined using a string containing C code), shell_volume must also be 
     966    #: C code, either defined as a string, or in the sources. 
     967    shell_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
    960968    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined 
    961969    #: by the parameter table.  *Iq* can be defined as a python function, or 
  • sasmodels/models/hollow_cylinder.c

    ree60aa7 re44432d  
    22 
    33static double 
    4 _fq(double qab, double qc, 
    5     double radius, double thickness, double length) 
     4shell_volume(double radius, double thickness, double length) 
    65{ 
    7     const double lam1 = sas_2J1x_x((radius+thickness)*qab); 
    8     const double lam2 = sas_2J1x_x(radius*qab); 
    9     const double gamma_sq = square(radius/(radius+thickness)); 
    10     //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab) 
    11     //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab) 
    12     const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00 
    13     const double t2 = sas_sinx_x(0.5*length*qc); 
    14     return psi*t2; 
    15 } 
    16  
    17 // TODO: interface to form_volume/shell_volume not yet settled 
    18 static double 
    19 shell_volume(double *total, double radius, double thickness, double length) 
    20 { 
    21     *total = M_PI*length*square(radius+thickness); 
    22     return *total - M_PI*length*radius*radius; 
     6    return M_PI*length*(square(radius+thickness) - radius*radius); 
    237} 
    248 
     
    2610form_volume(double radius, double thickness, double length) 
    2711{ 
    28     double total; 
    29     return shell_volume(&total, radius, thickness, length); 
     12    return M_PI*length*square(radius+thickness); 
    3013} 
    3114 
     
    6245} 
    6346 
     47static double 
     48_fq(double qab, double qc, 
     49    double radius, double thickness, double length) 
     50{ 
     51    const double lam1 = sas_2J1x_x((radius+thickness)*qab); 
     52    const double lam2 = sas_2J1x_x(radius*qab); 
     53    const double gamma_sq = square(radius/(radius+thickness)); 
     54    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab) 
     55    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab) 
     56    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00 
     57    const double t2 = sas_sinx_x(0.5*length*qc); 
     58    return psi*t2; 
     59} 
     60 
    6461static void 
    6562Fq(double q, double *F1, double *F2, double radius, double thickness, double length, 
     
    8178    total_F1 *= 0.5*(upper-lower); 
    8279    total_F2 *= 0.5*(upper-lower); 
    83     const double s = (sld - solvent_sld) * form_volume(radius, thickness, length); 
     80    const double s = (sld - solvent_sld) * shell_volume(radius, thickness, length); 
    8481    *F1 = 1e-2 * s * total_F1; 
    8582    *F2 = 1e-4 * s*s * total_F2; 
     
    9390{ 
    9491    const double form = _fq(qab, qc, radius, thickness, length); 
    95     const double s = (sld - solvent_sld) * form_volume(radius, thickness, length); 
     92    const double s = (sld - solvent_sld) * shell_volume(radius, thickness, length); 
    9693    return 1.0e-4*square(s * form); 
    9794} 
  • sasmodels/models/hollow_rectangular_prism.c

    ree60aa7 re44432d  
    1 // TODO: interface to form_volume/shell_volume not yet settled 
    21static double 
    3 shell_volume(double *total, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     2shell_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    43{ 
    5     double length_b = length_a * b2a_ratio; 
    6     double length_c = length_a * c2a_ratio; 
    7     double a_core = length_a - 2.0*thickness; 
    8     double b_core = length_b - 2.0*thickness; 
    9     double c_core = length_c - 2.0*thickness; 
    10     double vol_core = a_core * b_core * c_core; 
    11     *total = length_a * length_b * length_c; 
    12     return *total - vol_core; 
     4    const double length_b = length_a * b2a_ratio; 
     5    const double length_c = length_a * c2a_ratio; 
     6    const double form_volume = length_a * length_b * length_c; 
     7    const double a_core = length_a - 2.0*thickness; 
     8    const double b_core = length_b - 2.0*thickness; 
     9    const double c_core = length_c - 2.0*thickness; 
     10    const double core_volume = a_core * b_core * c_core; 
     11    return form_volume - core_volume; 
    1312} 
    1413 
     
    1615form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    1716{ 
    18     double total; 
    19     return shell_volume(&total, length_a, b2a_ratio, c2a_ratio, thickness); 
     17    const double length_b = length_a * b2a_ratio; 
     18    const double length_c = length_a * c2a_ratio; 
     19    const double form_volume = length_a * length_b * length_c; 
     20    return form_volume; 
    2021} 
    2122 
    2223static double 
    2324effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    24 //effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
    25 //                         "equivalent outer circular cross-section","half ab diagonal","half diagonal"] 
    2625// NOTE length_a is external dimension! 
    2726{ 
  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    ree60aa7 re44432d  
    1 // TODO: interface to form_volume/shell_volume not yet settled 
    21static double 
    3 shell_volume(double *total, double length_a, double b2a_ratio, double c2a_ratio) 
     2shell_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    43{ 
    5     double length_b = length_a * b2a_ratio; 
    6     double length_c = length_a * c2a_ratio; 
    7     double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
    8     *total = length_a * length_b * length_c; 
    9     return vol_shell; 
     4    const double length_b = length_a * b2a_ratio; 
     5    const double length_c = length_a * c2a_ratio; 
     6    const double shell_volume = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
     7    return shell_volume; 
    108} 
    119 
     
    1311form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    1412{ 
    15     double total; 
    16     return shell_volume(&total, length_a, b2a_ratio, c2a_ratio); 
     13    const double length_b = length_a * b2a_ratio; 
     14    const double length_c = length_a * c2a_ratio; 
     15    const double form_volume = length_a * length_b * length_c; 
     16    return form_volume; 
    1717} 
    18  
    1918 
    2019static double 
  • sasmodels/models/vesicle.c

    ree60aa7 re44432d  
    1 // TODO: interface to form_volume/shell_volume not yet settled 
    21static double 
    3 shell_volume(double *total, double radius, double thickness) 
     2shell_volume(double radius, double thickness) 
    43{ 
    5     //note that for the vesicle model, the volume is ONLY the shell volume 
    6     *total = M_4PI_3 * cube(radius+thickness); 
    7     return *total - M_4PI_3*cube(radius); 
     4    return M_4PI_3 * (cube(radius+thickness) - cube(radius)); 
    85} 
    96 
     
    118form_volume(double radius, double thickness) 
    129{ 
    13     //note that for the vesicle model, the volume is ONLY the shell volume 
    14     double total; 
    15     return shell_volume(&total, radius, thickness); 
     10    return M_4PI_3 * cube(radius+thickness); 
    1611} 
     12 
    1713 
    1814static double 
  • sasmodels/product.py

    r2773c66 re44432d  
    168168    return par 
    169169 
    170 def _intermediates(P, S, effective_radius): 
    171     # type: (np.ndarray, np.ndarray, float) -> OrderedDict[str, np.ndarray] 
    172     """ 
    173     Returns intermediate results for standard product (P(Q)*S(Q)) 
    174     """ 
    175     return OrderedDict(( 
    176         ("P(Q)", P), 
    177         ("S(Q)", S), 
    178         ("effective_radius", effective_radius), 
    179     )) 
    180  
    181 def _intermediates_beta(F1,              # type: np.ndarray 
    182                         F2,              # type: np.ndarray 
    183                         S,               # type: np.ndarray 
    184                         scale,           # type: np.ndarray 
    185                         bg,              # type: np.ndarray 
    186                         effective_radius # type: float 
    187                         ): 
     170def _intermediates( 
     171        F1,               # type: np.ndarray 
     172        F2,               # type: np.ndarray 
     173        S,                # type: np.ndarray 
     174        scale,            # type: float 
     175        effective_radius, # type: float 
     176        beta_mode,        # type: bool 
     177        ): 
    188178    # type: (...) -> OrderedDict[str, Union[np.ndarray, float]] 
    189179    """ 
    190     Returns intermediate results for beta approximation-enabled product. The result may be an array or a float. 
    191     """ 
    192     # TODO: 1. include calculated Q vector 
    193     # TODO: 2. consider implications if there are intermediate results in P(Q) 
    194     return OrderedDict(( 
    195         ("P(Q)", scale*F2), 
    196         ("S(Q)", S), 
    197         ("beta(Q)", F1**2 / F2), 
    198         ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), 
    199         ("effective_radius", effective_radius), 
    200         # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), 
    201     )) 
     180    Returns intermediate results for beta approximation-enabled product. 
     181    The result may be an array or a float. 
     182    """ 
     183    if beta_mode: 
     184        # TODO: 1. include calculated Q vector 
     185        # TODO: 2. consider implications if there are intermediate results in P(Q) 
     186        parts = OrderedDict(( 
     187            ("P(Q)", scale*F2), 
     188            ("S(Q)", S), 
     189            ("beta(Q)", F1**2 / F2), 
     190            ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), 
     191            ("effective_radius", effective_radius), 
     192            # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), 
     193        )) 
     194    else: 
     195        parts = OrderedDict(( 
     196            ("P(Q)", scale*F2), 
     197            ("S(Q)", S), 
     198            ("effective_radius", effective_radius), 
     199        )) 
     200    return parts 
    202201 
    203202class ProductModel(KernelModel): 
     
    253252    def __call__(self, call_details, values, cutoff, magnetic): 
    254253        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray 
     254 
    255255        p_info, s_info = self.info.composition[1] 
     256        p_npars = p_info.parameters.npars 
     257        p_length = call_details.length[:p_npars] 
     258        p_offset = call_details.offset[:p_npars] 
     259        s_npars = s_info.parameters.npars 
     260        s_length = call_details.length[p_npars:p_npars+s_npars] 
     261        s_offset = call_details.offset[p_npars:p_npars+s_npars] 
     262 
     263        # Beta mode parameter is the first parameter after P and S parameters 
     264        have_beta_mode = p_info.have_Fq 
     265        beta_mode_offset = 2+p_npars+s_npars 
     266        beta_mode = (values[beta_mode_offset] > 0) if have_beta_mode else False 
     267        if beta_mode and self.p_kernel.dim== '2d': 
     268            raise NotImplementedError("beta not yet supported for 2D") 
     269 
     270        # R_eff type parameter is the second parameter after P and S parameters 
     271        # unless the model doesn't support beta mode, in which case it is first 
     272        have_radius_type = p_info.effective_radius_type is not None 
     273        radius_type_offset = 2+p_npars+s_npars + (1 if have_beta_mode else 0) 
     274        radius_type = int(values[radius_type_offset]) if have_radius_type else 0 
     275 
     276        # Retrieve the volume fraction, which is the second of the 
     277        # 'S' parameters in the parameter list, or 2+np in 0-origin, 
     278        # as well as the scale and background. 
     279        volfrac = values[3+p_npars] 
     280        scale, background = values[0], values[1] 
    256281 
    257282        # if there are magnetic parameters, they will only be on the 
     
    268293 
    269294        # Construct the calling parameters for P. 
    270         p_npars = p_info.parameters.npars 
    271         p_length = call_details.length[:p_npars] 
    272         p_offset = call_details.offset[:p_npars] 
    273295        p_details = make_details(p_info, p_length, p_offset, nweights) 
    274  
    275         # Set p scale to the volume fraction in s, which is the first of the 
    276         # 'S' parameters in the parameter list, or 2+np in 0-origin. 
    277         volfrac = values[2+p_npars] 
    278         p_values = [[volfrac, 0.0], values[2:2+p_npars], magnetism, weights] 
     296        p_values = [ 
     297            [1., 0.], # scale=1, background=0, 
     298            values[2:2+p_npars], 
     299            magnetism, 
     300            weights] 
    279301        spacer = (32 - sum(len(v) for v in p_values)%32)%32 
    280302        p_values.append([0.]*spacer) 
     
    282304 
    283305        # Construct the calling parameters for S. 
    284         s_npars = s_info.parameters.npars 
    285         s_length = call_details.length[p_npars:p_npars+s_npars] 
    286         s_offset = call_details.offset[p_npars:p_npars+s_npars] 
    287         s_length = np.hstack((1, s_length)) 
    288         s_offset = np.hstack((nweights, s_offset)) 
    289         s_details = make_details(s_info, s_length, s_offset, nweights+1) 
     306        if radius_type > 0: 
     307            # If R_eff comes from form factor, make sure it is monodisperse. 
     308            # weight is set to 1 later, after the value array is created 
     309            s_length[0] = 1 
     310        s_details = make_details(s_info, s_length, s_offset, nweights) 
    290311        s_values = [ 
    291             # scale=1, background=0, 
    292             [1., 0.], 
     312            [1., 0.], # scale=1, background=0, 
    293313            values[2+p_npars:2+p_npars+s_npars], 
    294314            weights, 
     
    298318        s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 
    299319 
    300         # beta mode is the first parameter after the structure factor pars 
    301         extra_offset = 2+p_npars+s_npars 
    302         if p_info.have_Fq: 
    303             beta_mode = values[extra_offset] 
    304             extra_offset += 1 
    305         else: 
    306             beta_mode = 0 
    307         if p_info.effective_radius_type is not None: 
    308             effective_radius_type = int(values[extra_offset]) 
    309             extra_offset += 1 
    310         else: 
    311             effective_radius_type = 0 
    312  
    313         # Call the kernels 
    314         scale, background = values[0], values[1] 
    315         if beta_mode: 
    316             F1, F2, volume_avg, effective_radius = self.p_kernel.beta( 
    317                 p_details, p_values, cutoff, magnetic, effective_radius_type) 
    318             if effective_radius_type > 0: 
    319                 # Plug R_eff from p into S model (initial value and pd value) 
    320                 s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 
    321             s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    322             combined_scale = scale*volfrac/volume_avg 
    323  
    324             self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background, effective_radius) 
    325             final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background 
    326  
    327         else: 
    328             p_result, effective_radius = self.p_kernel.Pq_Reff( 
    329                 p_details, p_values, cutoff, magnetic, effective_radius_type) 
    330             if effective_radius_type > 0: 
    331                 # Plug R_eff from p into S model (initial value and pd value) 
    332                 s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 
    333             s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    334             # remember the parts for plotting later 
    335             self.results = lambda: _intermediates(p_result, s_result, effective_radius) 
    336             final_result = scale*(p_result*s_result) + background 
    337  
    338         #call_details.show(values) 
    339         #print("values", values) 
    340         #p_details.show(p_values) 
    341         #print("=>", p_result) 
    342         #s_details.show(s_values) 
    343         #print("=>", s_result) 
    344         #import pylab as plt 
    345         #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') 
    346         #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') 
    347         #plt.figure() 
     320        # Call the form factor kernel to compute <F> and <F^2>. 
     321        # If the model doesn't support Fq the returned <F> will be None. 
     322        F1, F2, effective_radius, shell_volume, volume_ratio = self.p_kernel.Fq( 
     323            p_details, p_values, cutoff, magnetic, radius_type) 
     324 
     325        # Call the structure factor kernel to compute S. 
     326        # Plug R_eff from the form factor into structure factor parameters 
     327        # and scale volume fraction by form:shell volume ratio. These changes 
     328        # needs to be both in the initial value slot as well as the 
     329        # polydispersity distribution slot in the values array due to 
     330        # implementation details in kernel_iq.c. 
     331        #print("R_eff=%d:%g, volfrac=%g, volume ratio=%g"%(radius_type, effective_radius, volfrac, volume_ratio)) 
     332        if radius_type > 0: 
     333            # set the value to the model ER and set the weight to 1 
     334            s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 
     335            s_values[2+s_npars+s_offset[0]+nweights] = 1.0 
     336        s_values[3] = s_values[2+s_npars+s_offset[1]] = volfrac*volume_ratio 
     337        S = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
     338 
     339        # Determine overall scale factor. Hollow shapes are weighted by 
     340        # shell_volume, so that is needed for volume normalization.  For 
     341        # solid shapes we can use shell_volume as well since it is equal 
     342        # to form volume. 
     343        combined_scale = scale*volfrac/shell_volume 
     344 
     345        # Combine form factor and structure factor 
     346        #print("beta", beta_mode, F1, F2, S) 
     347        PS = F2 + F1**2*(S-1) if beta_mode else F2*S 
     348        final_result = combined_scale*PS + background 
     349 
     350        # Capture intermediate values so user can see them.  These are 
     351        # returned as a lazy evaluator since they are only needed in the 
     352        # GUI, and not for each evaluation during a fit. 
     353        # TODO: return the results structure with the final results 
     354        # That way the model calcs are idempotent. Further, we can 
     355        # generalize intermediates to various other model types if we put it 
     356        # kernel calling interface.  Could do this as an "optional" 
     357        # return value in the caller, though in that case we could return 
     358        # the results directly rather than through a lazy evaluator. 
     359        self.results = lambda: _intermediates( 
     360            F1, F2, S, combined_scale, effective_radius, beta_mode) 
     361 
    348362        return final_result 
    349363 
  • sasmodels/sasview_model.py

    r9ae0d2e re44432d  
    847847    P = _make_standard_model('cylinder')() 
    848848    model = MultiplicationModel(P, S) 
     849    model.setParam('radius_effective_mode', 1.0) 
    849850    value = model.evalDistribution([0.1, 0.1]) 
    850851    if np.isnan(value): 
Note: See TracChangeset for help on using the changeset viewer.