Changeset e44432d in sasmodels
- Timestamp:
- Oct 19, 2018 5:46:26 PM (6 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- be43e39
- Parents:
- 2586ab72
- Location:
- sasmodels
- Files:
-
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/generate.py
r2773c66 re44432d 715 715 return False 716 716 717 _SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE) 718 def 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 717 729 def _add_source(source, code, path, lineno=1): 718 730 """ … … 767 779 pars = partable.form_volume_parameters 768 780 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)) 769 784 if isinstance(model_info.Iq, str): 770 785 pars = [q] + partable.iq_parameters … … 779 794 pars = [qa, qb, qc] + partable.iq_parameters 780 795 source.append(_gen_fn(model_info, 'Iqabc', pars)) 796 797 # Check for shell_volume in source 798 is_hollow = contains_shell_volume(source) 781 799 782 800 # What kind of 2D model do we need? Is it consistent with the parameters? … … 802 820 for p in partable.kernel_parameters)) 803 821 # 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" 805 823 if partable.form_volume_parameters: 806 824 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)) 808 829 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)) 810 831 else: 811 832 # Model doesn't have volume. We could make the kernel run a little 812 833 # faster by not using/transferring the volume normalizations, but 813 834 # 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)" 815 836 source.append(call_volume) 816 837 source.append(call_effective_radius) -
sasmodels/kernel.py
ree60aa7 re44432d 43 43 def Iq(self, call_details, values, cutoff, magnetic): 44 44 # 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 47 62 __call__ = Iq 48 63 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): 50 65 # 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 """ 51 131 self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 52 132 #print("returned",self.q_input.q, self.result) … … 55 135 if total_weight == 0.: 56 136 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 87 146 88 147 def release(self): … … 92 151 def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 93 152 # 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 """ 94 157 raise NotImplementedError() -
sasmodels/kernel_iq.c
r2773c66 re44432d 26 26 // MAGNETIC_PARS : a comma-separated list of indices to the sld 27 27 // parameters in the parameter table. 28 // CALL_VOLUME( table) : call the form volume function28 // CALL_VOLUME(form, shell, table) : assign form and shell values 29 29 // CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function 30 30 // CALL_IQ(q, table) : call the Iq function for 1D calcs. … … 275 275 276 276 // ==================== KERNEL CODE ======================== 277 #define COMPUTE_F1_F2 defined(CALL_FQ)278 277 kernel 279 278 void KERNEL_NAME( … … 345 344 // version must loop over all q. 346 345 #ifdef USE_OPENCL 347 #if COMPUTE_F1_F2346 #if defined(CALL_FQ) 348 347 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]); 351 351 double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]); 352 352 double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]); 353 353 #else 354 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]); 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]); 357 358 double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 358 359 #endif 359 360 #else // !USE_OPENCL 360 #if COMPUTE_F1_F2361 #if defined(CALL_FQ) 361 362 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]); 364 366 #else 365 367 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]); 368 371 #endif 369 372 if (pd_start == 0) { … … 371 374 #pragma omp parallel for 372 375 #endif 373 #if COMPUTE_F1_F2376 #if defined(CALL_FQ) 374 377 // 2*nq for F^2,F pairs 375 378 for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0; … … 378 381 #endif 379 382 } 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]); 381 384 #endif // !USE_OPENCL 382 385 … … 693 696 // Note: weight==0 must always be excluded 694 697 if (weight > cutoff) { 695 weighted_volume += weight * CALL_VOLUME(local_values.table);696 #if COMPUTE_F1_F2698 double form, shell; 699 CALL_VOLUME(form, shell, local_values.table); 697 700 weight_norm += weight; 698 #endif 701 weighted_form += weight * form; 702 weighted_shell += weight * shell; 699 703 if (effective_radius_type != 0) { 700 704 weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); … … 747 751 } 748 752 #else // !MAGNETIC 749 #if COMPUTE_F1_F2753 #if defined(CALL_FQ) 750 754 CALL_KERNEL(); // sets F1 and F2 by reference 751 755 #else … … 756 760 757 761 #ifdef USE_OPENCL 758 #if COMPUTE_F1_F2762 #if defined(CALL_FQ) 759 763 this_F2 += weight * F2; 760 764 this_F1 += weight * F1; … … 763 767 #endif 764 768 #else // !USE_OPENCL 765 #if COMPUTE_F1_F2769 #if defined(CALL_FQ) 766 770 result[2*q_index+0] += weight * F2; 767 771 result[2*q_index+1] += weight * F1; … … 793 797 // Remember the current result and the updated norm. 794 798 #ifdef USE_OPENCL 795 #if COMPUTE_F1_F2799 #if defined(CALL_FQ) 796 800 result[2*q_index+0] = this_F2; 797 801 result[2*q_index+1] = this_F1; 798 802 if (q_index == 0) { 799 803 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; 802 807 } 803 808 #else … … 805 810 if (q_index == 0) { 806 811 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; 809 815 } 810 816 #endif 811 817 812 //if (q_index == 0) printf("res: %g/%g\n", result[0], weig thed_volume);818 //if (q_index == 0) printf("res: %g/%g\n", result[0], weighted_shell); 813 819 #else // !USE_OPENCL 814 #if COMPUTE_F1_F2820 #if defined(CALL_FQ) 815 821 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; 818 825 #else 819 826 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; 822 830 #endif 823 //printf("res: %g/%g\n", result[0], weighted_ volume);831 //printf("res: %g/%g\n", result[0], weighted_shell); 824 832 #endif // !USE_OPENCL 825 833 826 834 // ** clear the macros in preparation for the next kernel ** 827 #undef COMPUTE_F1_F2828 835 #undef PD_INIT 829 836 #undef PD_OPEN -
sasmodels/kernelcl.py
r5399809 re44432d 539 539 # leave room for f1/f2 results in case we need to compute beta for 1d models 540 540 nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 541 # plus 3 weight, volume, radius542 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) 543 543 544 544 # Inputs and outputs for each kernel call -
sasmodels/kerneldll.py
r6e7ba14 re44432d 99 99 pass 100 100 # pylint: enable=unused-import 101 102 # Compiler output is a byte stream that needs to be decode in python 3 103 decode = (lambda s: s) if sys.version_info[0] < 3 else (lambda s: s.decode('utf8')) 101 104 102 105 if "SAS_DLL_PATH" in os.environ: … … 184 187 subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT) 185 188 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)) 187 191 if not os.path.exists(output): 188 192 raise RuntimeError("compile failed. File is in %r"%source) … … 383 387 # leave room for f1/f2 results in case we need to compute beta for 1d models 384 388 nout = 2 if self.info.have_Fq else 1 385 # plus 3 weight, volume, radius386 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) 387 391 self.real = (np.float32 if self.q_input.dtype == generate.F32 388 392 else np.float64 if self.q_input.dtype == generate.F64 -
sasmodels/kernelpy.py
r9ae0d2e re44432d 164 164 self._volume_args = volume_args 165 165 volume = model_info.form_volume 166 shell = model_info.shell_volume 166 167 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))) 169 171 self._radius = ((lambda mode: radius(mode, *volume_args)) if radius 170 172 else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume … … 225 227 total = form() 226 228 weight_norm = 1.0 227 weighted_ volume= form_volume()229 weighted_shell, weighted_form = form_volume() 228 230 weighted_radius = form_radius() 229 231 … … 233 235 234 236 weight_norm = 0.0 235 weighted_volume = 0.0 237 weighted_form = 0.0 238 weighted_shell = 0.0 236 239 weighted_radius = 0.0 237 240 partial_weight = np.NaN … … 272 275 total += weight * Iq 273 276 weight_norm += weight 274 weighted_volume += weight * form_volume() 277 shell, form = form_volume() 278 weighted_form += weight * form 279 weighted_shell += weight * shell 275 280 weighted_radius += weight * form_radius() 276 281 277 result = np.hstack((total, weight_norm, weighted_ volume, weighted_radius))282 result = np.hstack((total, weight_norm, weighted_form, weighted_shell, weighted_radius)) 278 283 return result 279 284 -
sasmodels/modelinfo.py
r9ae0d2e re44432d 718 718 719 719 #: Set of variables defined in the model that might contain C code 720 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', ' c_code']720 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code'] 721 721 722 722 def _find_source_lines(model_info, kernel_module): … … 803 803 info.tests = getattr(kernel_module, 'tests', []) 804 804 info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 805 info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore 805 806 info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 806 807 info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore … … 921 922 #: void Fq(double q, double *F1, double *F2, ...) 922 923 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] 923 927 #: List of C source files used to define the model. The source files 924 928 #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 925 929 #: model defines orientation parameters. Files containing the most basic 926 930 #: 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. 932 932 source = None # type: List[str] 933 933 #: The set of tests that must pass. The format of the tests is described … … 958 958 #: C code, either defined as a string, or in the sources. 959 959 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]] 960 968 #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined 961 969 #: by the parameter table. *Iq* can be defined as a python function, or -
sasmodels/models/hollow_cylinder.c
ree60aa7 re44432d 2 2 3 3 static double 4 _fq(double qab, double qc, 5 double radius, double thickness, double length) 4 shell_volume(double radius, double thickness, double length) 6 5 { 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); 23 7 } 24 8 … … 26 10 form_volume(double radius, double thickness, double length) 27 11 { 28 double total; 29 return shell_volume(&total, radius, thickness, length); 12 return M_PI*length*square(radius+thickness); 30 13 } 31 14 … … 62 45 } 63 46 47 static 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 64 61 static void 65 62 Fq(double q, double *F1, double *F2, double radius, double thickness, double length, … … 81 78 total_F1 *= 0.5*(upper-lower); 82 79 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); 84 81 *F1 = 1e-2 * s * total_F1; 85 82 *F2 = 1e-4 * s*s * total_F2; … … 93 90 { 94 91 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); 96 93 return 1.0e-4*square(s * form); 97 94 } -
sasmodels/models/hollow_rectangular_prism.c
ree60aa7 re44432d 1 // TODO: interface to form_volume/shell_volume not yet settled2 1 static double 3 shell_volume(double *total, doublelength_a, double b2a_ratio, double c2a_ratio, double thickness)2 shell_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 4 3 { 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; 13 12 } 14 13 … … 16 15 form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 17 16 { 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; 20 21 } 21 22 22 23 static double 23 24 effective_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"]26 25 // NOTE length_a is external dimension! 27 26 { -
sasmodels/models/hollow_rectangular_prism_thin_walls.c
ree60aa7 re44432d 1 // TODO: interface to form_volume/shell_volume not yet settled2 1 static double 3 shell_volume(double *total, doublelength_a, double b2a_ratio, double c2a_ratio)2 shell_volume(double length_a, double b2a_ratio, double c2a_ratio) 4 3 { 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; 10 8 } 11 9 … … 13 11 form_volume(double length_a, double b2a_ratio, double c2a_ratio) 14 12 { 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; 17 17 } 18 19 18 20 19 static double -
sasmodels/models/vesicle.c
ree60aa7 re44432d 1 // TODO: interface to form_volume/shell_volume not yet settled2 1 static double 3 shell_volume(double *total, doubleradius, double thickness)2 shell_volume(double radius, double thickness) 4 3 { 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)); 8 5 } 9 6 … … 11 8 form_volume(double radius, double thickness) 12 9 { 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); 16 11 } 12 17 13 18 14 static double -
sasmodels/product.py
r2773c66 re44432d 168 168 return par 169 169 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 ): 170 def _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 ): 188 178 # type: (...) -> OrderedDict[str, Union[np.ndarray, float]] 189 179 """ 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 202 201 203 202 class ProductModel(KernelModel): … … 253 252 def __call__(self, call_details, values, cutoff, magnetic): 254 253 # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray 254 255 255 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] 256 281 257 282 # if there are magnetic parameters, they will only be on the … … 268 293 269 294 # Construct the calling parameters for P. 270 p_npars = p_info.parameters.npars271 p_length = call_details.length[:p_npars]272 p_offset = call_details.offset[:p_npars]273 295 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 the276 # '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] 279 301 spacer = (32 - sum(len(v) for v in p_values)%32)%32 280 302 p_values.append([0.]*spacer) … … 282 304 283 305 # 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) 290 311 s_values = [ 291 # scale=1, background=0, 292 [1., 0.], 312 [1., 0.], # scale=1, background=0, 293 313 values[2+p_npars:2+p_npars+s_npars], 294 314 weights, … … 298 318 s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 299 319 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 348 362 return final_result 349 363 -
sasmodels/sasview_model.py
r9ae0d2e re44432d 847 847 P = _make_standard_model('cylinder')() 848 848 model = MultiplicationModel(P, S) 849 model.setParam('radius_effective_mode', 1.0) 849 850 value = model.evalDistribution([0.1, 0.1]) 850 851 if np.isnan(value):
Note: See TracChangeset
for help on using the changeset viewer.