Changes in / [62cf915:0ce5710] in sasmodels
- Files:
-
- 7 added
- 49 edited
Legend:
- Unmodified
- Added
- Removed
-
doc/developer/index.rst
r56fc97a rb85be2d 3 3 *************************** 4 4 5 .. toctree:: 6 :numbered: 4 7 :maxdepth: 4 5 8 9 calculator.rst -
doc/genapi.py
r3d5c6f8 ra5b8477 59 59 #('alignment', 'GPU data alignment [unused]'), 60 60 ('bumps_model', 'Bumps interface'), 61 ('compare', 'Compare models on different compute engines'), 61 62 ('convert', 'Sasview to sasmodel converter'), 62 63 ('core', 'Model access'), 64 ('data', 'Data layout and plotting routines'), 65 ('details', 'Parameter packing for kernel calls'), 63 66 ('direct_model', 'Simple interface'), 64 67 ('exception', 'Annotate exceptions'), 68 #('frozendict', 'Freeze a dictionary to make it immutable'), 65 69 ('generate', 'Model parser'), 70 ('kernel', 'Evaluator type definitions'), 66 71 ('kernelcl', 'OpenCL model evaluator'), 67 72 ('kerneldll', 'Ctypes model evaluator'), 68 73 ('kernelpy', 'Python model evaluator'), 74 ('list_pars', 'Identify all parameters in all models'), 75 ('mixture', 'Mixture model evaluator'), 69 76 ('model_test', 'Unit test support'), 77 ('modelinfo', 'Parameter and model definitions'), 78 ('product', 'Product model evaluator'), 70 79 ('resolution', '1-D resolution functions'), 71 80 ('resolution2d', '2-D resolution functions'), 72 81 ('sasview_model', 'Sasview interface'), 73 ('sesans', 'SESANS model evaluator'), 82 ('sesans', 'SESANS calculation routines'), 83 #('transition', 'Model stepper for automatic model selection'), 74 84 ('weights', 'Distribution functions'), 75 #('transition', 'Model stepper for automatic model selection'),76 85 ] 77 86 package='sasmodels' -
doc/genmodel.py
r89f4163 ra5b8477 1 from __future__ import print_function 2 1 3 import sys, os, math, re 2 4 import numpy as np 3 5 import matplotlib.pyplot as plt 4 import pylab5 6 sys.path.insert(0, os.path.abspath('..')) 6 7 from sasmodels import generate, core … … 8 9 from sasmodels.data import empty_data1D, empty_data2D 9 10 10 11 # Convert ../sasmodels/models/name.py to name 12 model_name = os.path.basename(sys.argv[1])[:-3] 13 model_info = core.load_model_info(model_name) 14 model = core.build_model(model_info) 15 16 # Load the doc string from the module definition file and store it in rst 17 docstr = generate.make_doc(model_info) 18 19 20 # Calculate 1D curve for default parameters 21 pars = dict((p.name, p.default) for p in model_info['parameters']) 22 23 # Plotting ranges and options 24 opts = { 25 'xscale' : 'log', 26 'yscale' : 'log' if not model_info['structure_factor'] else 'linear', 27 'zscale' : 'log' if not model_info['structure_factor'] else 'linear', 28 'q_min' : 0.001, 29 'q_max' : 1.0, 30 'nq' : 1000, 31 'nq2d' : 1000, 32 'vmin' : 1e-3, # floor for the 2D data results 33 'qx_max' : 0.5, 34 #'colormap' : 'gist_ncar', 35 'colormap' : 'nipy_spectral', 36 #'colormap' : 'jet', 37 } 38 11 try: 12 from typing import Dict, Any 13 except ImportError: 14 pass 15 else: 16 from matplotlib.axes import Axes 17 from sasmodels.kernel import KernelModel 18 from sasmodels.modelinfo import ModelInfo 39 19 40 20 def plot_1d(model, opts, ax): 21 # type: (KernelModel, Dict[str, Any], Axes) -> None 22 """ 23 Create a 1-D image. 24 """ 41 25 q_min, q_max, nq = opts['q_min'], opts['q_max'], opts['nq'] 42 26 q_min = math.log10(q_min) … … 47 31 Iq1D = calculator() 48 32 49 ax.plot(q, Iq1D, color='blue', lw=2, label=model _info['name'])33 ax.plot(q, Iq1D, color='blue', lw=2, label=model.info.name) 50 34 ax.set_xlabel(r'$Q \/(\AA^{-1})$') 51 35 ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$') … … 55 39 56 40 def plot_2d(model, opts, ax): 41 # type: (KernelModel, Dict[str, Any], Axes) -> None 42 """ 43 Create a 2-D image. 44 """ 57 45 qx_max, nq2d = opts['qx_max'], opts['nq2d'] 58 q = np.linspace(-qx_max, qx_max, nq2d) 46 q = np.linspace(-qx_max, qx_max, nq2d) # type: np.ndarray 59 47 data2d = empty_data2D(q, resolution=0.0) 60 48 calculator = DirectModel(data2d, model) … … 64 52 Iq2D = np.log(np.clip(Iq2D, opts['vmin'], np.inf)) 65 53 ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='lower', 66 extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=opts['colormap'])54 extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=opts['colormap']) 67 55 ax.set_xlabel(r'$Q_x \/(\AA^{-1})$') 68 56 ax.set_ylabel(r'$Q_y \/(\AA^{-1})$') 69 57 70 # Generate image 71 fig_height = 3.0 # in 72 fig_left = 0.6 # in 73 fig_right = 0.5 # in 74 fig_top = 0.6*0.25 # in 75 fig_bottom = 0.6*0.75 76 if model_info['has_2d']: 77 plot_height = fig_height - (fig_top+fig_bottom) 78 plot_width = plot_height 79 fig_width = 2*(plot_width + fig_left + fig_right) 80 aspect = (fig_width, fig_height) 81 ratio = aspect[0]/aspect[1] 82 ax_left = fig_left/fig_width 83 ax_bottom = fig_bottom/fig_height 84 ax_height = plot_height/fig_height 85 ax_width = ax_height/ratio # square axes 86 fig = plt.figure(figsize=aspect) 87 ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height]) 88 plot_2d(model, opts, ax2d) 89 ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 90 plot_1d(model, opts, ax1d) 91 #ax.set_aspect('square') 92 else: 93 plot_height = fig_height - (fig_top+fig_bottom) 94 plot_width = (1+np.sqrt(5))/2*fig_height 95 fig_width = plot_width + fig_left + fig_right 96 ax_left = fig_left/fig_width 97 ax_bottom = fig_bottom/fig_height 98 ax_width = plot_width/fig_width 99 ax_height = plot_height/fig_height 100 aspect = (fig_width, fig_height) 101 fig = plt.figure(figsize=aspect) 102 ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 103 plot_1d(model, opts, ax1d) 58 def figfile(model_info): 59 # type: (ModelInfo) -> str 60 return model_info.id + '_autogenfig.png' 104 61 105 # Save image in model/img 106 figname = model_name + '_autogenfig.png' 107 filename = os.path.join('model', 'img', figname) 108 plt.savefig(filename, bbox_inches='tight') 109 #print "figure saved in",filename 62 def make_figure(model_info, opts): 63 # type: (ModelInfo, Dict[str, Any]) -> None 64 """ 65 Generate the figure file to include in the docs. 66 """ 67 model = core.build_model(model_info) 110 68 111 # Auto caption for figure 112 captionstr = '\n' 113 captionstr += '.. figure:: img/' + model_info['id'] + '_autogenfig.png\n' 114 captionstr += '\n' 115 if model_info['has_2d']: 116 captionstr += ' 1D and 2D plots corresponding to the default parameters of the model.\n' 117 else: 118 captionstr += ' 1D plot corresponding to the default parameters of the model.\n' 119 captionstr += '\n' 69 fig_height = 3.0 # in 70 fig_left = 0.6 # in 71 fig_right = 0.5 # in 72 fig_top = 0.6*0.25 # in 73 fig_bottom = 0.6*0.75 74 if model_info.parameters.has_2d: 75 plot_height = fig_height - (fig_top+fig_bottom) 76 plot_width = plot_height 77 fig_width = 2*(plot_width + fig_left + fig_right) 78 aspect = (fig_width, fig_height) 79 ratio = aspect[0]/aspect[1] 80 ax_left = fig_left/fig_width 81 ax_bottom = fig_bottom/fig_height 82 ax_height = plot_height/fig_height 83 ax_width = ax_height/ratio # square axes 84 fig = plt.figure(figsize=aspect) 85 ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height]) 86 plot_2d(model, opts, ax2d) 87 ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 88 plot_1d(model, opts, ax1d) 89 #ax.set_aspect('square') 90 else: 91 plot_height = fig_height - (fig_top+fig_bottom) 92 plot_width = (1+np.sqrt(5))/2*fig_height 93 fig_width = plot_width + fig_left + fig_right 94 ax_left = fig_left/fig_width 95 ax_bottom = fig_bottom/fig_height 96 ax_width = plot_width/fig_width 97 ax_height = plot_height/fig_height 98 aspect = (fig_width, fig_height) 99 fig = plt.figure(figsize=aspect) 100 ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 101 plot_1d(model, opts, ax1d) 120 102 121 # Add figure reference and caption to documentation (at end, before References) 122 pattern = '\*\*REFERENCE' 123 m = re.search(pattern, docstr.upper()) 103 # Save image in model/img 104 path = os.path.join('model', 'img', figfile(model_info)) 105 plt.savefig(path, bbox_inches='tight') 106 #print("figure saved in",path) 124 107 125 if m: 126 docstr1 = docstr[:m.start()] 127 docstr2 = docstr[m.start():] 128 docstr = docstr1 + captionstr + docstr2 129 else: 130 print '------------------------------------------------------------------' 131 print 'References NOT FOUND for model: ', model_info['id'] 132 print '------------------------------------------------------------------' 133 docstr = docstr + captionstr 108 def gen_docs(model_info): 109 # type: (ModelInfo) -> None 110 """ 111 Generate the doc string with the figure inserted before the references. 112 """ 134 113 135 open(sys.argv[2],'w').write(docstr) 114 # Load the doc string from the module definition file and store it in rst 115 docstr = generate.make_doc(model_info) 116 117 # Auto caption for figure 118 captionstr = '\n' 119 captionstr += '.. figure:: img/' + figfile(model_info) + '\n' 120 captionstr += '\n' 121 if model_info.parameters.has_2d: 122 captionstr += ' 1D and 2D plots corresponding to the default parameters of the model.\n' 123 else: 124 captionstr += ' 1D plot corresponding to the default parameters of the model.\n' 125 captionstr += '\n' 126 127 # Add figure reference and caption to documentation (at end, before References) 128 pattern = '\*\*REFERENCE' 129 match = re.search(pattern, docstr.upper()) 130 131 if match: 132 docstr1 = docstr[:match.start()] 133 docstr2 = docstr[match.start():] 134 docstr = docstr1 + captionstr + docstr2 135 else: 136 print('------------------------------------------------------------------') 137 print('References NOT FOUND for model: ', model_info.id) 138 print('------------------------------------------------------------------') 139 docstr += captionstr 140 141 open(sys.argv[2],'w').write(docstr) 142 143 def process_model(path): 144 # type: (str) -> None 145 """ 146 Generate doc file and image file for the given model definition file. 147 """ 148 149 # Load the model file 150 model_name = os.path.basename(path)[:-3] 151 model_info = core.load_model_info(model_name) 152 153 # Plotting ranges and options 154 opts = { 155 'xscale' : 'log', 156 'yscale' : 'log' if not model_info.structure_factor else 'linear', 157 'zscale' : 'log' if not model_info.structure_factor else 'linear', 158 'q_min' : 0.001, 159 'q_max' : 1.0, 160 'nq' : 1000, 161 'nq2d' : 1000, 162 'vmin' : 1e-3, # floor for the 2D data results 163 'qx_max' : 0.5, 164 #'colormap' : 'gist_ncar', 165 'colormap' : 'nipy_spectral', 166 #'colormap' : 'jet', 167 } 168 169 # Generate the RST file and the figure. Order doesn't matter. 170 gen_docs(model_info) 171 make_figure(model_info, opts) 172 173 if __name__ == "__main__": 174 process_model(sys.argv[1]) -
doc/gentoc.py
r5041682 ra5b8477 9 9 from sasmodels.core import load_model_info 10 10 11 try: 12 from typing import Optional, BinaryIO, List, Dict 13 except ImportError: 14 pass 15 else: 16 from sasmodels.modelinfo import ModelInfo 11 17 12 18 TEMPLATE="""\ … … 27 33 28 34 def _make_category(category_name, label, title, parent=None): 35 # type: (str, str, str, Optional[BinaryIO]) -> BinaryIO 29 36 file = open(joinpath(MODEL_TOC_PATH, category_name+".rst"), "w") 30 37 file.write(TEMPLATE%{'label':label, 'title':title, 'bar':'*'*len(title)}) … … 34 41 35 42 def _add_subcategory(category_name, parent): 43 # type: (str, BinaryIO) -> None 36 44 parent.write(" %s.rst\n"%category_name) 37 45 38 46 def _add_model(file, model_name): 47 # type: (IO[str], str) -> None 39 48 file.write(" ../../model/%s.rst\n"%model_name) 40 49 41 50 def _maybe_make_category(category, models, cat_files, model_toc): 51 # type: (str, List[str], Dict[str, BinaryIO], BinaryIO) -> None 42 52 if category not in cat_files: 43 53 print("Unexpected category %s containing"%category, models, file=sys.stderr) … … 46 56 47 57 def generate_toc(model_files): 58 # type: (List[str]) -> None 48 59 if not model_files: 49 60 print("gentoc needs a list of model files", file=sys.stderr) 50 61 51 62 # find all categories 52 category = {} 63 category = {} # type: Dict[str, List[str]] 53 64 for item in model_files: 54 65 # assume model is in sasmodels/models/name.py, and ignore the full path … … 56 67 if model_name.startswith('_'): continue 57 68 model_info = load_model_info(model_name) 58 if model_info ['category']is None:69 if model_info.category is None: 59 70 print("Missing category for", item, file=sys.stderr) 60 71 else: 61 category.setdefault(model_info ['category'],[]).append(model_name)72 category.setdefault(model_info.category,[]).append(model_name) 62 73 63 74 # Check category names -
sasmodels/alignment.py
r3c56da87 r7ae2b7f 18 18 to decide whether it is really required. 19 19 """ 20 import numpy as np 20 import numpy as np # type: ignore 21 21 22 22 def align_empty(shape, dtype, alignment=128): -
sasmodels/bumps_model.py
rea75043 r04dc697 11 11 12 12 """ 13 14 import warnings 15 16 import numpy as np 13 from __future__ import print_function 14 15 __all__ = [ "Model", "Experiment" ] 16 17 import numpy as np # type: ignore 17 18 18 19 from .data import plot_theory 19 20 from .direct_model import DataMixin 20 21 21 __all__ = [ 22 "Model", "Experiment", 23 ] 24 25 # CRUFT: old style bumps wrapper which doesn't separate data and model 26 # pylint: disable=invalid-name 27 def BumpsModel(data, model, cutoff=1e-5, **kw): 28 r""" 29 Bind a model to data, along with a polydispersity cutoff. 30 31 *data* is a :class:`data.Data1D`, :class:`data.Data2D` or 32 :class:`data.Sesans` object. Use :func:`data.empty_data1D` or 33 :func:`data.empty_data2D` to define $q, \Delta q$ calculation 34 points for displaying the SANS curve when there is no measured data. 35 36 *model* is a runnable module as returned from :func:`core.load_model`. 37 38 *cutoff* is the polydispersity weight cutoff. 39 40 Any additional *key=value* pairs are model dependent parameters. 41 42 Returns an :class:`Experiment` object. 43 44 Note that the usual Bumps semantics is not fully supported, since 45 assigning *M.name = parameter* on the returned experiment object 46 does not set that parameter in the model. Range setting will still 47 work as expected though. 48 49 .. deprecated:: 0.1 50 Use :class:`Experiment` instead. 51 """ 52 warnings.warn("Use of BumpsModel is deprecated. Use bumps_model.Experiment instead.") 53 54 # Create the model and experiment 55 model = Model(model, **kw) 56 experiment = Experiment(data=data, model=model, cutoff=cutoff) 57 58 # Copy the model parameters up to the experiment object. 59 for k, v in model.parameters().items(): 60 setattr(experiment, k, v) 61 return experiment 22 try: 23 from typing import Dict, Union, Tuple, Any 24 from .data import Data1D, Data2D 25 from .kernel import KernelModel 26 from .modelinfo import ModelInfo 27 Data = Union[Data1D, Data2D] 28 except ImportError: 29 pass 30 31 try: 32 # Optional import. This allows the doc builder and nosetests to run even 33 # when bumps is not on the path. 34 from bumps.names import Parameter # type: ignore 35 except ImportError: 36 pass 62 37 63 38 64 39 def create_parameters(model_info, **kwargs): 40 # type: (ModelInfo, **Union[float, str, Parameter]) -> Tuple[Dict[str, Parameter], Dict[str, str]] 65 41 """ 66 42 Generate Bumps parameters from the model info. … … 71 47 Any additional *key=value* pairs are initial values for the parameters 72 48 to the models. Uninitialized parameters will use the model default 73 value. 49 value. The value can be a float, a bumps parameter, or in the case 50 of the distribution type parameter, a string. 74 51 75 52 Returns a dictionary of *{name: Parameter}* containing the bumps … … 77 54 *{name: str}* containing the polydispersity distribution types. 78 55 """ 79 # lazy import; this allows the doc builder and nosetests to run even 80 # when bumps is not on the path. 81 from bumps.names import Parameter 82 83 pars = {} 84 for p in model_info['parameters']: 56 pars = {} # type: Dict[str, Parameter] 57 pd_types = {} # type: Dict[str, str] 58 for p in model_info.parameters.call_parameters: 85 59 value = kwargs.pop(p.name, p.default) 86 60 pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits) 87 for name in model_info['partype']['pd-2d']: 88 for xpart, xdefault, xlimits in [ 89 ('_pd', 0., pars[name].limits), 90 ('_pd_n', 35., (0, 1000)), 91 ('_pd_nsigma', 3., (0, 10)), 92 ]: 93 xname = name + xpart 94 xvalue = kwargs.pop(xname, xdefault) 95 pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) 96 97 pd_types = {} 98 for name in model_info['partype']['pd-2d']: 99 xname = name + '_pd_type' 100 xvalue = kwargs.pop(xname, 'gaussian') 101 pd_types[xname] = xvalue 61 if p.polydisperse: 62 for part, default, limits in [ 63 ('_pd', 0., pars[p.name].limits), 64 ('_pd_n', 35., (0, 1000)), 65 ('_pd_nsigma', 3., (0, 10)), 66 ]: 67 name = p.name + part 68 value = kwargs.pop(name, default) 69 pars[name] = Parameter.default(value, name=name, limits=limits) 70 name = p.name + '_pd_type' 71 pd_types[name] = str(kwargs.pop(name, 'gaussian')) 102 72 103 73 if kwargs: # args not corresponding to parameters … … 118 88 """ 119 89 def __init__(self, model, **kwargs): 120 self._sasmodel = model 90 # type: (KernelModel, **Dict[str, Union[float, Parameter]]) -> None 91 self.sasmodel = model 121 92 pars, pd_types = create_parameters(model.info, **kwargs) 122 93 for k, v in pars.items(): … … 128 99 129 100 def parameters(self): 101 # type: () -> Dict[str, Parameter] 130 102 """ 131 103 Return a dictionary of parameters objects for the parameters, … … 135 107 136 108 def state(self): 109 # type: () -> Dict[str, Union[Parameter, str]] 137 110 """ 138 111 Return a dictionary of current values for all the parameters, … … 159 132 The resulting model can be used directly in a Bumps FitProblem call. 160 133 """ 134 _cache = None # type: Dict[str, np.ndarray] 161 135 def __init__(self, data, model, cutoff=1e-5): 162 136 # type: (Data, Model, float) -> None 163 137 # remember inputs so we can inspect from outside 164 138 self.model = model 165 139 self.cutoff = cutoff 166 self._interpret_data(data, model. _sasmodel)167 self. update()140 self._interpret_data(data, model.sasmodel) 141 self._cache = {} 168 142 169 143 def update(self): 144 # type: () -> None 170 145 """ 171 146 Call when model parameters have changed and theory needs to be 172 147 recalculated. 173 148 """ 174 self._cache = {}149 self._cache.clear() 175 150 176 151 def numpoints(self): 152 # type: () -> float 177 153 """ 178 154 Return the number of data points … … 181 157 182 158 def parameters(self): 159 # type: () -> Dict[str, Parameter] 183 160 """ 184 161 Return a dictionary of parameters … … 187 164 188 165 def theory(self): 166 # type: () -> np.ndarray 189 167 """ 190 168 Return the theory corresponding to the model parameters. … … 199 177 200 178 def residuals(self): 179 # type: () -> np.ndarray 201 180 """ 202 181 Return theory minus data normalized by uncertainty. … … 206 185 207 186 def nllf(self): 187 # type: () -> float 208 188 """ 209 189 Return the negative log likelihood of seeing data given the model … … 213 193 delta = self.residuals() 214 194 #if np.any(np.isnan(R)): print("NaN in residuals") 215 return 0.5 * np.sum(delta **2)195 return 0.5 * np.sum(delta**2) 216 196 217 197 #def __call__(self): … … 219 199 220 200 def plot(self, view='log'): 201 # type: (str) -> None 221 202 """ 222 203 Plot the data and residuals. … … 226 207 227 208 def simulate_data(self, noise=None): 209 # type: (float) -> None 228 210 """ 229 211 Generate simulated data. … … 233 215 234 216 def save(self, basename): 217 # type: (str) -> None 235 218 """ 236 219 Save the model parameters and data into a file. … … 243 226 244 227 def __getstate__(self): 228 # type: () -> Dict[str, Any] 245 229 # Can't pickle gpu functions, so instead make them lazy 246 230 state = self.__dict__.copy() … … 249 233 250 234 def __setstate__(self, state): 235 # type: (Dict[str, Any]) -> None 251 236 # pylint: disable=attribute-defined-outside-init 252 237 self.__dict__ = state 238 -
sasmodels/compare.py
rf247314 r7ae2b7f 34 34 import traceback 35 35 36 import numpy as np 36 import numpy as np # type: ignore 37 37 38 38 from . import core 39 39 from . import kerneldll 40 from . import product41 40 from .data import plot_theory, empty_data1D, empty_data2D 42 41 from .direct_model import DirectModel … … 210 209 Choose a parameter range based on parameter name and initial value. 211 210 """ 211 # process the polydispersity options 212 212 if p.endswith('_pd_n'): 213 213 return [0, 100] … … 222 222 else: 223 223 return [-180, 180] 224 elif p.endswith('_pd'): 225 return [0, 1] 224 226 elif 'sld' in p: 225 227 return [-0.5, 10] 226 elif p.endswith('_pd'):227 return [0, 1]228 228 elif p == 'background': 229 229 return [0, 10] 230 230 elif p == 'scale': 231 231 return [0, 1e3] 232 elif p == 'case_num':233 # RPA hack234 return [0, 10]235 232 elif v < 0: 236 # Kxy parameters in rpa model can be negative237 233 return [2*v, -2*v] 238 234 else: … … 240 236 241 237 242 def _randomize_one( p, v):238 def _randomize_one(model_info, p, v): 243 239 """ 244 240 Randomize a single parameter. … … 246 242 if any(p.endswith(s) for s in ('_pd_n', '_pd_nsigma', '_pd_type')): 247 243 return v 244 245 # Find the parameter definition 246 for par in model_info.parameters.call_parameters: 247 if par.name == p: 248 break 248 249 else: 249 return np.random.uniform(*parameter_range(p, v)) 250 251 252 def randomize_pars(pars, seed=None): 250 raise ValueError("unknown parameter %r"%p) 251 252 if len(par.limits) > 2: # choice list 253 return np.random.randint(len(par.limits)) 254 255 limits = par.limits 256 if not np.isfinite(limits).all(): 257 low, high = parameter_range(p, v) 258 limits = (max(limits[0], low), min(limits[1], high)) 259 260 return np.random.uniform(*limits) 261 262 263 def randomize_pars(model_info, pars, seed=None): 253 264 """ 254 265 Generate random values for all of the parameters. … … 261 272 with push_seed(seed): 262 273 # Note: the sort guarantees order `of calls to random number generator 263 pars = dict((p, _randomize_one( p, v))274 pars = dict((p, _randomize_one(model_info, p, v)) 264 275 for p, v in sorted(pars.items())) 265 276 return pars … … 273 284 cylinder radius in this case). 274 285 """ 275 name = model_info ['id']286 name = model_info.id 276 287 # if it is a product model, then just look at the form factor since 277 288 # none of the structure factors need any constraints. … … 294 305 # Make sure phi sums to 1.0 295 306 if pars['case_num'] < 2: 296 pars['Phi a'] = 0.297 pars['Phi b'] = 0.307 pars['Phi1'] = 0. 308 pars['Phi2'] = 0. 298 309 elif pars['case_num'] < 5: 299 pars['Phi a'] = 0.300 total = sum(pars['Phi'+c] for c in ' abcd')301 for c in ' abcd':310 pars['Phi1'] = 0. 311 total = sum(pars['Phi'+c] for c in '1234') 312 for c in '1234': 302 313 pars['Phi'+c] /= total 303 314 … … 306 317 Format the parameter list for printing. 307 318 """ 308 if is2d:309 exclude = lambda n: False310 else:311 partype = model_info['partype']312 par1d = set(partype['fixed-1d']+partype['pd-1d'])313 exclude = lambda n: n not in par1d314 319 lines = [] 315 for p in model_info['parameters']:316 if exclude(p.name): continue320 parameters = model_info.parameters 321 for p in parameters.user_parameters(pars, is2d): 317 322 fields = dict( 318 value=pars.get(p. name, p.default),319 pd=pars.get(p. name+"_pd", 0.),320 n=int(pars.get(p. name+"_pd_n", 0)),321 nsigma=pars.get(p. name+"_pd_nsgima", 3.),322 type=pars.get(p. name+"_pd_type", 'gaussian'))323 value=pars.get(p.id, p.default), 324 pd=pars.get(p.id+"_pd", 0.), 325 n=int(pars.get(p.id+"_pd_n", 0)), 326 nsigma=pars.get(p.id+"_pd_nsgima", 3.), 327 type=pars.get(p.id+"_pd_type", 'gaussian')) 323 328 lines.append(_format_par(p.name, **fields)) 324 329 return "\n".join(lines) … … 363 368 364 369 # grab the sasview model, or create it if it is a product model 365 if model_info ['composition']:366 composition_type, parts = model_info ['composition']370 if model_info.composition: 371 composition_type, parts = model_info.composition 367 372 if composition_type == 'product': 368 373 from sas.models.MultiplicationModel import MultiplicationModel … … 454 459 """ 455 460 # initialize the code so time is more accurate 456 value = calculator(**suppress_pd(pars)) 461 if Nevals > 1: 462 value = calculator(**suppress_pd(pars)) 457 463 toc = tic() 458 464 for _ in range(max(Nevals, 1)): # make sure there is at least one eval … … 661 667 """ 662 668 # Get the default values for the parameters 663 pars = dict((p.name, p.default) for p in model_info['parameters']) 664 665 # Fill in default values for the polydispersity parameters 666 for p in model_info['parameters']: 667 if p.type in ('volume', 'orientation'): 668 pars[p.name+'_pd'] = 0.0 669 pars[p.name+'_pd_n'] = 0 670 pars[p.name+'_pd_nsigma'] = 3.0 671 pars[p.name+'_pd_type'] = "gaussian" 669 pars = {} 670 for p in model_info.parameters.call_parameters: 671 parts = [('', p.default)] 672 if p.polydisperse: 673 parts.append(('_pd', 0.0)) 674 parts.append(('_pd_n', 0)) 675 parts.append(('_pd_nsigma', 3.0)) 676 parts.append(('_pd_type', "gaussian")) 677 for ext, val in parts: 678 if p.length > 1: 679 dict(("%s%d%s"%(p.id,k,ext), val) for k in range(p.length)) 680 else: 681 pars[p.id+ext] = val 672 682 673 683 # Plug in values given in demo 674 684 if use_demo: 675 pars.update(model_info ['demo'])685 pars.update(model_info.demo) 676 686 return pars 677 687 … … 700 710 try: 701 711 model_info = core.load_model_info(name) 702 except ImportError ,exc:712 except ImportError as exc: 703 713 print(str(exc)) 704 714 print("Could not find model; use one of:\n " + models) … … 774 784 775 785 if len(engines) == 0: 776 engines.extend(['single', ' sasview'])786 engines.extend(['single', 'double']) 777 787 elif len(engines) == 1: 778 if engines[0][0] != ' sasview':779 engines.append(' sasview')788 if engines[0][0] != 'double': 789 engines.append('double') 780 790 else: 781 791 engines.append('single') … … 807 817 #pars.update(set_pars) # set value before random to control range 808 818 if opts['seed'] > -1: 809 pars = randomize_pars( pars, seed=opts['seed'])819 pars = randomize_pars(model_info, pars, seed=opts['seed']) 810 820 print("Randomize using -random=%i"%opts['seed']) 811 821 if opts['mono']: … … 850 860 Explore the model using the Bumps GUI. 851 861 """ 852 import wx 853 from bumps.names import FitProblem 854 from bumps.gui.app_frame import AppFrame 862 import wx # type: ignore 863 from bumps.names import FitProblem # type: ignore 864 from bumps.gui.app_frame import AppFrame # type: ignore 855 865 856 866 problem = FitProblem(Explore(opts)) … … 873 883 """ 874 884 def __init__(self, opts): 875 from bumps.cli import config_matplotlib 885 from bumps.cli import config_matplotlib # type: ignore 876 886 from . import bumps_model 877 887 config_matplotlib() … … 879 889 model_info = opts['def'] 880 890 pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 891 # Initialize parameter ranges, fixing the 2D parameters for 1D data. 881 892 if not opts['is2d']: 882 active = [base + ext 883 for base in model_info['partype']['pd-1d'] 884 for ext in ['', '_pd', '_pd_n', '_pd_nsigma']] 885 active.extend(model_info['partype']['fixed-1d']) 886 for k in active: 887 v = pars[k] 888 v.range(*parameter_range(k, v.value)) 893 for p in model_info.parameters.user_parameters(is2d=False): 894 for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 895 k = p.name+ext 896 v = pars.get(k, None) 897 if v is not None: 898 v.range(*parameter_range(k, v.value)) 889 899 else: 890 900 for k, v in pars.items(): -
sasmodels/compare_many.py
rce346b6 r7ae2b7f 17 17 import traceback 18 18 19 import numpy as np 19 import numpy as np # type: ignore 20 20 21 21 from . import core 22 from . import generate23 22 from .compare import (MODELS, randomize_pars, suppress_pd, make_data, 24 23 make_engine, get_pars, columnize, … … 109 108 110 109 if is_2d: 111 if not model_info[' has_2d']:110 if not model_info['parameters'].has_2d: 112 111 print(',"1-D only"') 113 112 return … … 125 124 try: 126 125 result = fn(**pars) 127 except KeyboardInterrupt: 128 raise 129 except: 126 except Exception: 130 127 traceback.print_exc() 131 128 print("when comparing %s for %d"%(name, seed)) … … 246 243 base = sys.argv[5] 247 244 comp = sys.argv[6] if len(sys.argv) > 6 else "sasview" 248 except :245 except Exception: 249 246 traceback.print_exc() 250 247 print_usage() -
sasmodels/convert.json
rec45c4f r6e7ff6d 602 602 "RPAModel", 603 603 { 604 "N1": "Na", "N2": "Nb", "N3": "Nc", "N4": "Nd", 605 "L1": "La", "L2": "Lb", "L3": "Lc", "L4": "Ld", 606 "v1": "va", "v2": "vb", "v3": "vc", "v4": "vd", 607 "b1": "ba", "b2": "bb", "b3": "bc", "b4": "bd", 608 "Phi1": "Phia", "Phi2": "Phib", "Phi3": "Phic", "Phi4": "Phid", 604 609 "case_num": "lcase_n" 605 610 } … … 620 625 } 621 626 ], 627 "_spherepy": [ 628 "SphereModel", 629 { 630 "sld": "sldSph", 631 "radius": "radius", 632 "sld_solvent": "sldSolv" 633 } 634 ], 622 635 "spherical_sld": [ 623 636 "SphereSLDModel", 624 637 { 638 "n": "n_shells", 625 639 "radius_core": "rad_core", 626 640 "sld_solvent": "sld_solv" -
sasmodels/convert.py
rf247314 r7ae2b7f 62 62 # but only if we haven't already byteified it 63 63 if isinstance(data, dict) and not ignore_dicts: 64 return { 65 _byteify(key, ignore_dicts=True): _byteify(value, ignore_dicts=True) 66 for key, value in data.iteritems() 67 } 64 return dict((_byteify(key, ignore_dicts=True), 65 _byteify(value, ignore_dicts=True)) 66 for key, value in data.items()) 68 67 # if it's anything else, return it in its original form 69 68 return data … … 153 152 def revert_name(model_info): 154 153 _read_conversion_table() 155 oldname, oldpars = CONVERSION_TABLE.get(model_info ['id'], [None, {}])154 oldname, oldpars = CONVERSION_TABLE.get(model_info.id, [None, {}]) 156 155 return oldname 157 156 158 157 def _get_old_pars(model_info): 159 158 _read_conversion_table() 160 oldname, oldpars = CONVERSION_TABLE.get(model_info ['id'], [None, {}])159 oldname, oldpars = CONVERSION_TABLE.get(model_info.id, [None, {}]) 161 160 return oldpars 162 161 … … 165 164 Convert model from new style parameter names to old style. 166 165 """ 167 if model_info ['composition']is not None:168 composition_type, parts = model_info ['composition']166 if model_info.composition is not None: 167 composition_type, parts = model_info.composition 169 168 if composition_type == 'product': 170 169 P, S = parts … … 180 179 181 180 # Note: update compare.constrain_pars to match 182 name = model_info ['id']183 if name in MODELS_WITHOUT_SCALE or model_info ['structure_factor']:181 name = model_info.id 182 if name in MODELS_WITHOUT_SCALE or model_info.structure_factor: 184 183 if oldpars.pop('scale', 1.0) != 1.0: 185 184 warnings.warn("parameter scale not used in sasview %s"%name) 186 if name in MODELS_WITHOUT_BACKGROUND or model_info ['structure_factor']:185 if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor: 187 186 if oldpars.pop('background', 0.0) != 0.0: 188 187 warnings.warn("parameter background not used in sasview %s"%name) … … 206 205 elif name == 'rpa': 207 206 # convert scattering lengths from femtometers to centimeters 208 for p in "L a", "Lb", "Lc", "Ld":207 for p in "L1", "L2", "L3", "L4": 209 208 if p in oldpars: oldpars[p] *= 1e-13 210 209 elif name == 'core_shell_parallelepiped': … … 225 224 Restrict parameter values to those that will match sasview. 226 225 """ 227 name = model_info ['id']226 name = model_info.id 228 227 # Note: update convert.revert_model to match 229 if name in MODELS_WITHOUT_SCALE or model_info ['structure_factor']:228 if name in MODELS_WITHOUT_SCALE or model_info.structure_factor: 230 229 pars['scale'] = 1 231 if name in MODELS_WITHOUT_BACKGROUND or model_info ['structure_factor']:230 if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor: 232 231 pars['background'] = 0 233 232 # sasview multiplies background by structure factor -
sasmodels/core.py
r4d76711 rfa5fd8d 2 2 Core model handling routines. 3 3 """ 4 from __future__ import print_function 5 4 6 __all__ = [ 5 "list_models", "load_model _info", "precompile_dll",6 "build_model", " make_kernel", "call_kernel", "call_ER_VR",7 "list_models", "load_model", "load_model_info", 8 "build_model", "precompile_dll", 7 9 ] 8 10 9 from os.path import basename, dirname, join as joinpath , splitext11 from os.path import basename, dirname, join as joinpath 10 12 from glob import glob 11 13 12 import numpy as np 14 import numpy as np # type: ignore 13 15 14 from . import models15 from . import weights16 16 from . import generate 17 # TODO: remove circular references between product and core 18 # product uses call_ER/call_VR, core uses make_product_info/ProductModel 19 #from . import product 17 from . import modelinfo 18 from . import product 20 19 from . import mixture 21 20 from . import kernelpy … … 24 23 from . import kernelcl 25 24 HAVE_OPENCL = True 26 except :25 except Exception: 27 26 HAVE_OPENCL = False 28 27 29 28 try: 30 np.meshgrid([]) 31 meshgrid = np.meshgrid 32 except ValueError: 33 # CRUFT: np.meshgrid requires multiple vectors 34 def meshgrid(*args): 35 if len(args) > 1: 36 return np.meshgrid(*args) 37 else: 38 return [np.asarray(v) for v in args] 29 from typing import List, Union, Optional, Any 30 DType = Union[None, str, np.dtype] 31 from .kernel import KernelModel 32 except ImportError: 33 pass 34 39 35 40 36 # TODO: refactor composite model support … … 51 47 52 48 def list_models(): 49 # type: () -> List[str] 53 50 """ 54 51 Return the list of available models on the model path. … … 60 57 61 58 def isstr(s): 59 # type: (Any) -> bool 62 60 """ 63 61 Return True if *s* is a string-like object. 64 62 """ 65 63 try: s + '' 66 except : return False64 except Exception: return False 67 65 return True 68 66 69 def load_model(model_name, **kw): 67 def load_model(model_name, dtype=None, platform='ocl'): 68 # type: (str, DType, str) -> KernelModel 70 69 """ 71 70 Load model info and build model. 71 72 *model_name* is the name of the model as used by :func:`load_model_info`. 73 Additional keyword arguments are passed directly to :func:`build_model`. 72 74 """ 73 return build_model(load_model_info(model_name), **kw) 75 return build_model(load_model_info(model_name), 76 dtype=dtype, platform=platform) 74 77 75 78 76 79 def load_model_info(model_name): 80 # type: (str) -> modelinfo.ModelInfo 77 81 """ 78 82 Load a model definition given the model name. … … 88 92 parts = model_name.split('*') 89 93 if len(parts) > 1: 90 from . import product91 # Note: currently have circular reference92 94 if len(parts) > 2: 93 95 raise ValueError("use P*S to apply structure factor S to model P") … … 96 98 97 99 kernel_module = generate.load_kernel_module(model_name) 98 return generate.make_model_info(kernel_module)100 return modelinfo.make_model_info(kernel_module) 99 101 100 102 101 103 def build_model(model_info, dtype=None, platform="ocl"): 104 # type: (modelinfo.ModelInfo, np.dtype, str) -> KernelModel 102 105 """ 103 106 Prepare the model for the default execution platform. … … 118 121 otherwise it uses the default "ocl". 119 122 """ 120 composition = model_info. get('composition', None)123 composition = model_info.composition 121 124 if composition is not None: 122 125 composition_type, parts = composition … … 131 134 raise ValueError('unknown mixture type %s'%composition_type) 132 135 136 # If it is a python model, return it immediately 137 if callable(model_info.Iq): 138 return kernelpy.PyModel(model_info) 139 133 140 ## for debugging: 134 141 ## 1. uncomment open().write so that the source will be saved next time … … 137 144 ## 4. rerun "python -m sasmodels.direct_model $MODELNAME" 138 145 ## 5. uncomment open().read() so that source will be regenerated from model 139 # open(model_info ['name']+'.c','w').write(source)140 # source = open(model_info ['name']+'.cl','r').read()146 # open(model_info.name+'.c','w').write(source) 147 # source = open(model_info.name+'.cl','r').read() 141 148 source = generate.make_source(model_info) 142 149 if dtype is None: 143 dtype = 'single' if model_info['single'] else 'double' 144 if callable(model_info.get('Iq', None)): 145 return kernelpy.PyModel(model_info) 150 dtype = generate.F32 if model_info.single else generate.F64 146 151 if (platform == "dll" 147 152 or not HAVE_OPENCL … … 152 157 153 158 def precompile_dll(model_name, dtype="double"): 159 # type: (str, DType) -> Optional[str] 154 160 """ 155 161 Precompile the dll for a model. … … 168 174 source = generate.make_source(model_info) 169 175 return kerneldll.make_dll(source, model_info, dtype=dtype) if source else None 170 171 172 def get_weights(model_info, pars, name):173 """174 Generate the distribution for parameter *name* given the parameter values175 in *pars*.176 177 Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma"178 from the *pars* dictionary for parameter value and parameter dispersion.179 """180 relative = name in model_info['partype']['pd-rel']181 limits = model_info['limits'][name]182 disperser = pars.get(name+'_pd_type', 'gaussian')183 value = pars.get(name, model_info['defaults'][name])184 npts = pars.get(name+'_pd_n', 0)185 width = pars.get(name+'_pd', 0.0)186 nsigma = pars.get(name+'_pd_nsigma', 3.0)187 value, weight = weights.get_weights(188 disperser, npts, width, nsigma, value, limits, relative)189 return value, weight / np.sum(weight)190 191 def dispersion_mesh(pars):192 """193 Create a mesh grid of dispersion parameters and weights.194 195 Returns [p1,p2,...],w where pj is a vector of values for parameter j196 and w is a vector containing the products for weights for each197 parameter set in the vector.198 """199 value, weight = zip(*pars)200 value = [v.flatten() for v in meshgrid(*value)]201 weight = np.vstack([v.flatten() for v in meshgrid(*weight)])202 weight = np.prod(weight, axis=0)203 return value, weight204 205 def call_kernel(kernel, pars, cutoff=0, mono=False):206 """207 Call *kernel* returned from *model.make_kernel* with parameters *pars*.208 209 *cutoff* is the limiting value for the product of dispersion weights used210 to perform the multidimensional dispersion calculation more quickly at a211 slight cost to accuracy. The default value of *cutoff=0* integrates over212 the entire dispersion cube. Using *cutoff=1e-5* can be 50% faster, but213 with an error of about 1%, which is usually less than the measurement214 uncertainty.215 216 *mono* is True if polydispersity should be set to none on all parameters.217 """218 fixed_pars = [pars.get(name, kernel.info['defaults'][name])219 for name in kernel.fixed_pars]220 if mono:221 pd_pars = [( np.array([pars[name]]), np.array([1.0]) )222 for name in kernel.pd_pars]223 else:224 pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars]225 return kernel(fixed_pars, pd_pars, cutoff=cutoff)226 227 def call_ER_VR(model_info, vol_pars):228 """229 Return effect radius and volume ratio for the model.230 231 *info* is either *kernel.info* for *kernel=make_kernel(model,q)*232 or *model.info*.233 234 *pars* are the parameters as expected by :func:`call_kernel`.235 """236 ER = model_info.get('ER', None)237 VR = model_info.get('VR', None)238 value, weight = dispersion_mesh(vol_pars)239 240 individual_radii = ER(*value) if ER else 1.0241 whole, part = VR(*value) if VR else (1.0, 1.0)242 243 effect_radius = np.sum(weight*individual_radii) / np.sum(weight)244 volume_ratio = np.sum(weight*part)/np.sum(weight*whole)245 return effect_radius, volume_ratio246 247 248 def call_ER(model_info, values):249 """250 Call the model ER function using *values*. *model_info* is either251 *model.info* if you have a loaded model, or *kernel.info* if you252 have a model kernel prepared for evaluation.253 """254 ER = model_info.get('ER', None)255 if ER is None:256 return 1.0257 else:258 vol_pars = [get_weights(model_info, values, name)259 for name in model_info['partype']['volume']]260 value, weight = dispersion_mesh(vol_pars)261 individual_radii = ER(*value)262 #print(values[0].shape, weights.shape, fv.shape)263 return np.sum(weight*individual_radii) / np.sum(weight)264 265 def call_VR(model_info, values):266 """267 Call the model VR function using *pars*.268 *info* is either *model.info* if you have a loaded model, or *kernel.info*269 if you have a model kernel prepared for evaluation.270 """271 VR = model_info.get('VR', None)272 if VR is None:273 return 1.0274 else:275 vol_pars = [get_weights(model_info, values, name)276 for name in model_info['partype']['volume']]277 value, weight = dispersion_mesh(vol_pars)278 whole, part = VR(*value)279 return np.sum(weight*part)/np.sum(weight*whole)280 281 # TODO: remove call_ER, call_VR282 -
sasmodels/custom/__init__.py
rcd0a808 r7ae2b7f 14 14 try: 15 15 # Python 3.5 and up 16 from importlib.util import spec_from_file_location, module_from_spec 16 from importlib.util import spec_from_file_location, module_from_spec # type: ignore 17 17 def load_module_from_path(fullname, path): 18 18 spec = spec_from_file_location(fullname, path) -
sasmodels/data.py
rd6f5da6 ra5b8477 35 35 import traceback 36 36 37 import numpy as np 37 import numpy as np # type: ignore 38 39 try: 40 from typing import Union, Dict, List, Optional 41 except ImportError: 42 pass 43 else: 44 Data = Union["Data1D", "Data2D", "SesansData"] 38 45 39 46 def load_data(filename): 47 # type: (str) -> Data 40 48 """ 41 49 Load data using a sasview loader. 42 50 """ 43 from sas.sascalc.dataloader.loader import Loader 51 from sas.sascalc.dataloader.loader import Loader # type: ignore 44 52 loader = Loader() 45 53 data = loader.load(filename) … … 50 58 51 59 def set_beam_stop(data, radius, outer=None): 60 # type: (Data, float, Optional[float]) -> None 52 61 """ 53 62 Add a beam stop of the given *radius*. If *outer*, make an annulus. 54 63 """ 55 from sas.dataloader.manipulations import Ringcut 64 from sas.dataloader.manipulations import Ringcut # type: ignore 56 65 if hasattr(data, 'qx_data'): 57 66 data.mask = Ringcut(0, radius)(data) … … 65 74 66 75 def set_half(data, half): 76 # type: (Data, str) -> None 67 77 """ 68 78 Select half of the data, either "right" or "left". 69 79 """ 70 from sas.dataloader.manipulations import Boxcut 80 from sas.dataloader.manipulations import Boxcut # type: ignore 71 81 if half == 'right': 72 82 data.mask += \ … … 78 88 79 89 def set_top(data, cutoff): 90 # type: (Data, float) -> None 80 91 """ 81 92 Chop the top off the data, above *cutoff*. 82 93 """ 83 from sas.dataloader.manipulations import Boxcut 94 from sas.dataloader.manipulations import Boxcut # type: ignore 84 95 data.mask += \ 85 96 Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=cutoff)(data) … … 114 125 """ 115 126 def __init__(self, x=None, y=None, dx=None, dy=None): 127 # type: (Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray]) -> None 116 128 self.x, self.y, self.dx, self.dy = x, y, dx, dy 117 129 self.dxl = None … … 127 139 128 140 def xaxis(self, label, unit): 141 # type: (str, str) -> None 129 142 """ 130 143 set the x axis label and unit … … 134 147 135 148 def yaxis(self, label, unit): 149 # type: (str, str) -> None 136 150 """ 137 151 set the y axis label and unit … … 140 154 self._yunit = unit 141 155 142 156 class SesansData(Data1D): 157 def __init__(self, **kw): 158 Data1D.__init__(self, **kw) 159 self.lam = None # type: Optional[np.ndarray] 143 160 144 161 class Data2D(object): … … 175 192 """ 176 193 def __init__(self, x=None, y=None, z=None, dx=None, dy=None, dz=None): 194 # type: (Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray]) -> None 177 195 self.qx_data, self.dqx_data = x, dx 178 196 self.qy_data, self.dqy_data = y, dy … … 197 215 198 216 def xaxis(self, label, unit): 217 # type: (str, str) -> None 199 218 """ 200 219 set the x axis label and unit … … 204 223 205 224 def yaxis(self, label, unit): 225 # type: (str, str) -> None 206 226 """ 207 227 set the y axis label and unit … … 211 231 212 232 def zaxis(self, label, unit): 233 # type: (str, str) -> None 213 234 """ 214 235 set the y axis label and unit … … 223 244 """ 224 245 def __init__(self, x=None, y=None, z=None): 246 # type: (float, float, Optional[float]) -> None 225 247 self.x, self.y, self.z = x, y, z 226 248 … … 230 252 """ 231 253 def __init__(self, pixel_size=(None, None), distance=None): 254 # type: (Tuple[float, float], float) -> None 232 255 self.pixel_size = Vector(*pixel_size) 233 256 self.distance = distance … … 238 261 """ 239 262 def __init__(self): 263 # type: () -> None 240 264 self.wavelength = np.NaN 241 265 self.wavelength_unit = "A" … … 243 267 244 268 def empty_data1D(q, resolution=0.0): 269 # type: (np.ndarray, float) -> Data1D 245 270 """ 246 271 Create empty 1D data using the given *q* as the x value. … … 259 284 260 285 def empty_data2D(qx, qy=None, resolution=0.0): 286 # type: (np.ndarray, Optional[np.ndarray], float) -> Data2D 261 287 """ 262 288 Create empty 2D data using the given mesh. … … 272 298 Qx, Qy = np.meshgrid(qx, qy) 273 299 Qx, Qy = Qx.flatten(), Qy.flatten() 274 Iq = 100 * np.ones_like(Qx) 300 Iq = 100 * np.ones_like(Qx) # type: np.ndarray 275 301 dIq = np.sqrt(Iq) 276 302 if resolution != 0: … … 300 326 301 327 def plot_data(data, view='log', limits=None): 328 # type: (Data, str, Optional[Tuple[float, float]]) -> None 302 329 """ 303 330 Plot data loaded by the sasview loader. … … 323 350 def plot_theory(data, theory, resid=None, view='log', 324 351 use_data=True, limits=None, Iq_calc=None): 352 # type: (Data, Optional[np.ndarray], Optional[np.ndarray], str, bool, Optional[Tuple[float,float]], Optional[np.ndarray]) -> None 325 353 """ 326 354 Plot theory calculation. … … 337 365 *limits* sets the intensity limits on the plot; if None then the limits 338 366 are inferred from the data. 367 368 *Iq_calc* is the raw theory values without resolution smearing 339 369 """ 340 370 if hasattr(data, 'lam'): … … 348 378 349 379 def protect(fn): 380 # type: (Callable) -> Callable 350 381 """ 351 382 Decorator to wrap calls in an exception trapper which prints the … … 358 389 try: 359 390 return fn(*args, **kw) 360 except KeyboardInterrupt: 361 raise 362 except: 391 except Exception: 363 392 traceback.print_exc() 364 393 … … 369 398 def _plot_result1D(data, theory, resid, view, use_data, 370 399 limits=None, Iq_calc=None): 400 # type: (Data1D, Optional[np.ndarray], Optional[np.ndarray], str, bool, Optional[Tuple[float, float]], Optional[np.ndarray]) -> None 371 401 """ 372 402 Plot the data and residuals for 1D data. 373 403 """ 374 import matplotlib.pyplot as plt 375 from numpy.ma import masked_array, masked 404 import matplotlib.pyplot as plt # type: ignore 405 from numpy.ma import masked_array, masked # type: ignore 376 406 377 407 use_data = use_data and data.y is not None … … 446 476 @protect 447 477 def _plot_result_sesans(data, theory, resid, use_data, limits=None): 478 # type: (SesansData, Optional[np.ndarray], Optional[np.ndarray], bool, Optional[Tuple[float, float]]) -> None 448 479 """ 449 480 Plot SESANS results. 450 481 """ 451 import matplotlib.pyplot as plt 482 import matplotlib.pyplot as plt # type: ignore 452 483 use_data = use_data and data.y is not None 453 484 use_theory = theory is not None … … 456 487 457 488 if use_data or use_theory: 458 is_tof = np.any(data.lam!=data.lam[0])489 is_tof = (data.lam != data.lam[0]).any() 459 490 if num_plots > 1: 460 491 plt.subplot(1, num_plots, 1) 461 492 if use_data: 462 493 if is_tof: 463 plt.errorbar(data.x, np.log(data.y)/(data.lam*data.lam), yerr=data.dy/data.y/(data.lam*data.lam)) 494 plt.errorbar(data.x, np.log(data.y)/(data.lam*data.lam), 495 yerr=data.dy/data.y/(data.lam*data.lam)) 464 496 else: 465 497 plt.errorbar(data.x, data.y, yerr=data.dy) … … 489 521 @protect 490 522 def _plot_result2D(data, theory, resid, view, use_data, limits=None): 523 # type: (Data2D, Optional[np.ndarray], Optional[np.ndarray], str, bool, Optional[Tuple[float,float]]) -> None 491 524 """ 492 525 Plot the data and residuals for 2D data. 493 526 """ 494 import matplotlib.pyplot as plt 527 import matplotlib.pyplot as plt # type: ignore 495 528 use_data = use_data and data.data is not None 496 529 use_theory = theory is not None … … 500 533 # Put theory and data on a common colormap scale 501 534 vmin, vmax = np.inf, -np.inf 535 target = None # type: Optional[np.ndarray] 502 536 if use_data: 503 537 target = data.data[~data.mask] … … 548 582 @protect 549 583 def _plot_2d_signal(data, signal, vmin=None, vmax=None, view='log'): 584 # type: (Data2D, np.ndarray, Optional[float], Optional[float], str) -> Tuple[float, float] 550 585 """ 551 586 Plot the target value for the data. This could be the data itself, … … 554 589 *scale* can be 'log' for log scale data, or 'linear'. 555 590 """ 556 import matplotlib.pyplot as plt 557 from numpy.ma import masked_array 591 import matplotlib.pyplot as plt # type: ignore 592 from numpy.ma import masked_array # type: ignore 558 593 559 594 image = np.zeros_like(data.qx_data) … … 589 624 590 625 def demo(): 626 # type: () -> None 591 627 """ 592 628 Load and plot a SAS dataset. … … 595 631 set_beam_stop(data, 0.004) 596 632 plot_data(data) 597 import matplotlib.pyplot as plt; plt.show() 633 import matplotlib.pyplot as plt # type: ignore 634 plt.show() 598 635 599 636 -
sasmodels/direct_model.py
r4d76711 ra5b8477 23 23 from __future__ import print_function 24 24 25 import numpy as np 26 27 from .core import call_kernel, call_ER_VR 28 from . import sesans 25 import numpy as np # type: ignore 26 27 # TODO: fix sesans module 28 from . import sesans # type: ignore 29 from . import weights 29 30 from . import resolution 30 31 from . import resolution2d 32 from . import details 33 34 try: 35 from typing import Optional, Dict, Tuple 36 except ImportError: 37 pass 38 else: 39 from .data import Data 40 from .kernel import Kernel, KernelModel 41 from .modelinfo import Parameter, ParameterSet 42 43 def call_kernel(kernel, pars, cutoff=0., mono=False): 44 # type: (Kernel, ParameterSet, float, bool) -> np.ndarray 45 """ 46 Call *kernel* returned from *model.make_kernel* with parameters *pars*. 47 48 *cutoff* is the limiting value for the product of dispersion weights used 49 to perform the multidimensional dispersion calculation more quickly at a 50 slight cost to accuracy. The default value of *cutoff=0* integrates over 51 the entire dispersion cube. Using *cutoff=1e-5* can be 50% faster, but 52 with an error of about 1%, which is usually less than the measurement 53 uncertainty. 54 55 *mono* is True if polydispersity should be set to none on all parameters. 56 """ 57 parameters = kernel.info.parameters 58 if mono: 59 active = lambda name: False 60 elif kernel.dim == '1d': 61 active = lambda name: name in parameters.pd_1d 62 elif kernel.dim == '2d': 63 active = lambda name: name in parameters.pd_2d 64 else: 65 active = lambda name: True 66 67 vw_pairs = [(get_weights(p, pars) if active(p.name) 68 else ([pars.get(p.name, p.default)], [])) 69 for p in parameters.call_parameters] 70 71 call_details, weights, values = details.build_details(kernel, vw_pairs) 72 return kernel(call_details, weights, values, cutoff) 73 74 def get_weights(parameter, values): 75 # type: (Parameter, Dict[str, float]) -> Tuple[np.ndarray, np.ndarray] 76 """ 77 Generate the distribution for parameter *name* given the parameter values 78 in *pars*. 79 80 Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma" 81 from the *pars* dictionary for parameter value and parameter dispersion. 82 """ 83 value = float(values.get(parameter.name, parameter.default)) 84 relative = parameter.relative_pd 85 limits = parameter.limits 86 disperser = values.get(parameter.name+'_pd_type', 'gaussian') 87 npts = values.get(parameter.name+'_pd_n', 0) 88 width = values.get(parameter.name+'_pd', 0.0) 89 nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 90 if npts == 0 or width == 0: 91 return [value], [] 92 value, weight = weights.get_weights( 93 disperser, npts, width, nsigma, value, limits, relative) 94 return value, weight / np.sum(weight) 31 95 32 96 class DataMixin(object): … … 52 116 """ 53 117 def _interpret_data(self, data, model): 118 # type: (Data, KernelModel) -> None 54 119 # pylint: disable=attribute-defined-outside-init 55 120 … … 67 132 self.data_type = 'Iq' 68 133 69 partype = model.info['partype']70 71 134 if self.data_type == 'sesans': 72 135 … … 82 145 q_mono = sesans.make_all_q(data) 83 146 elif self.data_type == 'Iqxy': 84 if not partype['orientation'] and not partype['magnetic']:85 raise ValueError("not 2D without orientation or magnetic parameters")147 #if not model.info.parameters.has_2d: 148 # raise ValueError("not 2D without orientation or magnetic parameters") 86 149 q = np.sqrt(data.qx_data**2 + data.qy_data**2) 87 150 qmin = getattr(data, 'qmin', 1e-16) … … 154 217 155 218 def _set_data(self, Iq, noise=None): 219 # type: (np.ndarray, Optional[float]) -> None 156 220 # pylint: disable=attribute-defined-outside-init 157 221 if noise is not None: … … 171 235 172 236 def _calc_theory(self, pars, cutoff=0.0): 237 # type: (ParameterSet, float) -> np.ndarray 173 238 if self._kernel is None: 174 self._kernel = self._model.make_kernel(self._kernel_inputs) # pylint: disable=attribute-dedata_type239 self._kernel = self._model.make_kernel(self._kernel_inputs) 175 240 self._kernel_mono = (self._model.make_kernel(self._kernel_mono_inputs) 176 241 if self._kernel_mono_inputs else None) … … 206 271 """ 207 272 def __init__(self, data, model, cutoff=1e-5): 273 # type: (Data, KernelModel, float) -> None 208 274 self.model = model 209 275 self.cutoff = cutoff … … 212 278 213 279 def __call__(self, **pars): 280 # type: (**float) -> np.ndarray 214 281 return self._calc_theory(pars, cutoff=self.cutoff) 215 282 216 def ER_VR(self, **pars):217 """218 Compute the equivalent radius and volume ratio for the model.219 """220 return call_ER_VR(self.model.info, pars)221 222 283 def simulate_data(self, noise=None, **pars): 284 # type: (Optional[float], **float) -> None 223 285 """ 224 286 Generate simulated data for the model. … … 228 290 229 291 def main(): 292 # type: () -> None 230 293 """ 231 294 Program to evaluate a particular model at a set of q values. … … 243 306 try: 244 307 values = [float(v) for v in call.split(',')] 245 except :308 except Exception: 246 309 values = [] 247 310 if len(values) == 1: -
sasmodels/generate.py
r81ec7c8 ra5b8477 21 21 22 22 *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms. 23 24 #define INVALID(v) (expr) returns False if v.parameter is invalid 25 for some parameter or other (e.g., v.bell_radius < v.radius). If 26 necessary, the expression can call a function. 23 27 24 28 These functions are defined in a kernel module .py script and an associated … … 65 69 The constructor code will generate the necessary vectors for computing 66 70 them with the desired polydispersity. 67 68 The available kernel parameters are defined as a list, with each parameter69 defined as a sublist with the following elements:70 71 *name* is the name that will be used in the call to the kernel72 function and the name that will be displayed to the user. Names73 should be lower case, with words separated by underscore. If74 acronyms are used, the whole acronym should be upper case.75 76 *units* should be one of *degrees* for angles, *Ang* for lengths,77 *1e-6/Ang^2* for SLDs.78 79 *default value* will be the initial value for the model when it80 is selected, or when an initial value is not otherwise specified.81 82 *limits = [lb, ub]* are the hard limits on the parameter value, used to83 limit the polydispersity density function. In the fit, the parameter limits84 given to the fit are the limits on the central value of the parameter.85 If there is polydispersity, it will evaluate parameter values outside86 the fit limits, but not outside the hard limits specified in the model.87 If there are no limits, use +/-inf imported from numpy.88 89 *type* indicates how the parameter will be used. "volume" parameters90 will be used in all functions. "orientation" parameters will be used91 in *Iqxy* and *Imagnetic*. "magnetic* parameters will be used in92 *Imagnetic* only. If *type* is the empty string, the parameter will93 be used in all of *Iq*, *Iqxy* and *Imagnetic*. "sld" parameters94 can automatically be promoted to magnetic parameters, each of which95 will have a magnitude and a direction, which may be different from96 other sld parameters.97 98 *description* is a short description of the parameter. This will99 be displayed in the parameter table and used as a tool tip for the100 parameter value in the user interface.101 102 71 The kernel module must set variables defining the kernel meta data: 103 72 … … 154 123 polydispersity values for tests. 155 124 156 An *model_info* dictionary is constructed from the kernel meta data and 157 returned to the caller. 158 159 The model evaluator, function call sequence consists of q inputs and the return vector, 160 followed by the loop value/weight vector, followed by the values for 161 the non-polydisperse parameters, followed by the lengths of the 162 polydispersity loops. To construct the call for 1D models, the 163 categories *fixed-1d* and *pd-1d* list the names of the parameters 164 of the non-polydisperse and the polydisperse parameters respectively. 165 Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models. 166 The *pd-rel* category is a set of those parameters which give 167 polydispersitiy as a portion of the value (so a 10% length dispersity 168 would use a polydispersity value of 0.1) rather than absolute 169 dispersity such as an angle plus or minus 15 degrees. 170 171 The *volume* category lists the volume parameters in order for calls 172 to volume within the kernel (used for volume normalization) and for 173 calls to ER and VR for effective radius and volume ratio respectively. 174 175 The *orientation* and *magnetic* categories list the orientation and 176 magnetic parameters. These are used by the sasview interface. The 177 blank category is for parameters such as scale which don't have any 178 other marking. 125 A :class:`modelinfo.ModelInfo` structure is constructed from the kernel meta 126 data and returned to the caller. 179 127 180 128 The doc string at the start of the kernel module will be used to … … 184 132 file systems are case-sensitive, so only use lower case characters for 185 133 file names and extensions. 186 187 188 The function :func:`make` loads the metadata from the module and returns189 the kernel source. The function :func:`make_doc` extracts the doc string190 and adds the parameter table to the top. The function :func:`model_sources`191 returns a list of files required by the model.192 134 193 135 Code follows the C99 standard with the following extensions and conditions:: … … 201 143 all double declarations may be converted to half, float, or long double 202 144 FLOAT_SIZE is the number of bytes in the converted variables 145 146 :func:`load_kernel_module` loads the model definition file and 147 :modelinfo:`make_model_info` parses it. :func:`make_source` 148 converts C-based model definitions to C source code, including the 149 polydispersity integral. :func:`model_sources` returns the list of 150 source files the model depends on, and :func:`timestamp` returns 151 the latest time stamp amongst the source files (so you can check if 152 the model needs to be rebuilt). 153 154 The function :func:`make_doc` extracts the doc string and adds the 155 parameter table to the top. *make_figure* in *sasmodels/doc/genmodel* 156 creates the default figure for the model. [These two sets of code 157 should mignrate into docs.py so docs can be updated in one place]. 203 158 """ 204 159 from __future__ import print_function 205 160 206 #TODO: identify model files which have changed since loading and reload them.207 161 #TODO: determine which functions are useful outside of generate 208 162 #__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 209 163 210 import sys 211 from os.path import abspath, dirname, join as joinpath, exists, basename, \ 212 splitext 164 from os.path import abspath, dirname, join as joinpath, exists, getmtime 213 165 import re 214 166 import string 215 167 import warnings 216 from collections import namedtuple 217 218 import numpy as np 219 168 169 import numpy as np # type: ignore 170 171 from .modelinfo import Parameter 220 172 from .custom import load_custom_kernel_module 221 173 222 PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description'] 223 Parameter = namedtuple('Parameter', PARAMETER_FIELDS) 224 225 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 174 try: 175 from typing import Tuple, Sequence, Iterator 176 from .modelinfo import ModelInfo 177 except ImportError: 178 pass 179 180 TEMPLATE_ROOT = dirname(__file__) 226 181 227 182 F16 = np.dtype('float16') … … 232 187 except TypeError: 233 188 F128 = None 234 235 # Scale and background, which are parameters common to every form factor236 COMMON_PARAMETERS = [237 ["scale", "", 1, [0, np.inf], "", "Source intensity"],238 ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"],239 ]240 189 241 190 # Conversion from units defined in the parameter table for each model … … 284 233 285 234 def format_units(units): 235 # type: (str) -> str 286 236 """ 287 237 Convert units into ReStructured Text format. … … 290 240 291 241 def make_partable(pars): 242 # type: (List[Parameter]) -> str 292 243 """ 293 244 Generate the parameter table to include in the sphinx documentation. … … 320 271 321 272 def _search(search_path, filename): 273 # type: (List[str], str) -> str 322 274 """ 323 275 Find *filename* in *search_path*. … … 331 283 raise ValueError("%r not found in %s" % (filename, search_path)) 332 284 285 333 286 def model_sources(model_info): 287 # type: (ModelInfo) -> List[str] 334 288 """ 335 289 Return a list of the sources file paths for the module. 336 290 """ 337 search_path = [dirname(model_info ['filename']),291 search_path = [dirname(model_info.filename), 338 292 abspath(joinpath(dirname(__file__), 'models'))] 339 return [_search(search_path, f) for f in model_info['source']] 340 341 # Pragmas for enable OpenCL features. Be sure to protect them so that they 342 # still compile even if OpenCL is not present. 343 _F16_PRAGMA = """\ 344 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 345 # pragma OPENCL EXTENSION cl_khr_fp16: enable 346 #endif 347 """ 348 349 _F64_PRAGMA = """\ 350 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 351 # pragma OPENCL EXTENSION cl_khr_fp64: enable 352 #endif 353 """ 293 return [_search(search_path, f) for f in model_info.source] 294 295 def timestamp(model_info): 296 # type: (ModelInfo) -> int 297 """ 298 Return a timestamp for the model corresponding to the most recently 299 changed file or dependency. 300 """ 301 source_files = (model_sources(model_info) 302 + model_templates() 303 + [model_info.filename]) 304 newest = max(getmtime(f) for f in source_files) 305 return newest 354 306 355 307 def convert_type(source, dtype): 308 # type: (str, np.dtype) -> str 356 309 """ 357 310 Convert code from double precision to the desired type. … … 362 315 if dtype == F16: 363 316 fbytes = 2 364 source = _ F16_PRAGMA + _convert_type(source, "half", "f")317 source = _convert_type(source, "half", "f") 365 318 elif dtype == F32: 366 319 fbytes = 4 … … 368 321 elif dtype == F64: 369 322 fbytes = 8 370 source = _F64_PRAGMA + source # Sourceis already double323 # no need to convert if it is already double 371 324 elif dtype == F128: 372 325 fbytes = 16 … … 378 331 379 332 def _convert_type(source, type_name, constant_flag): 333 # type: (str, str, str) -> str 380 334 """ 381 335 Replace 'double' with *type_name* in *source*, tagging floating point … … 396 350 397 351 def kernel_name(model_info, is_2d): 352 # type: (ModelInfo, bool) -> str 398 353 """ 399 354 Name of the exported kernel symbol. 400 355 """ 401 return model_info ['name']+ "_" + ("Iqxy" if is_2d else "Iq")356 return model_info.name + "_" + ("Iqxy" if is_2d else "Iq") 402 357 403 358 404 359 def indent(s, depth): 360 # type: (str, int) -> str 405 361 """ 406 362 Indent a string of text with *depth* additional spaces on each line. … … 411 367 412 368 413 LOOP_OPEN = """\ 414 for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 415 const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 416 const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 369 _template_cache = {} # type: Dict[str, Tuple[int, str, str]] 370 def load_template(filename): 371 # type: (str) -> str 372 path = joinpath(TEMPLATE_ROOT, filename) 373 mtime = getmtime(path) 374 if filename not in _template_cache or mtime > _template_cache[filename][0]: 375 with open(path) as fid: 376 _template_cache[filename] = (mtime, fid.read(), path) 377 return _template_cache[filename][1] 378 379 def model_templates(): 380 # type: () -> List[str] 381 # TODO: fails DRY; templates are listed in two places. 382 # should instead have model_info contain a list of paths 383 return [joinpath(TEMPLATE_ROOT, filename) 384 for filename in ('kernel_header.c', 'kernel_iq.c')] 385 386 387 _FN_TEMPLATE = """\ 388 double %(name)s(%(pars)s); 389 double %(name)s(%(pars)s) { 390 %(body)s 391 } 392 393 417 394 """ 418 def build_polydispersity_loops(pd_pars): 419 """ 420 Build polydispersity loops 421 422 Returns loop opening and loop closing 423 """ 424 depth = 4 425 offset = "" 426 loop_head = [] 427 loop_end = [] 428 for name in pd_pars: 429 subst = {'name': name, 'offset': offset} 430 loop_head.append(indent(LOOP_OPEN % subst, depth)) 431 loop_end.insert(0, (" "*depth) + "}") 432 offset += '+N' + name 433 depth += 2 434 return "\n".join(loop_head), "\n".join(loop_end) 435 436 C_KERNEL_TEMPLATE = None 395 396 def _gen_fn(name, pars, body): 397 # type: (str, List[Parameter], str) -> str 398 """ 399 Generate a function given pars and body. 400 401 Returns the following string:: 402 403 double fn(double a, double b, ...); 404 double fn(double a, double b, ...) { 405 .... 406 } 407 """ 408 par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 409 return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 410 411 def _call_pars(prefix, pars): 412 # type: (str, List[Parameter]) -> List[str] 413 """ 414 Return a list of *prefix.parameter* from parameter items. 415 """ 416 return [p.as_call_reference(prefix) for p in pars] 417 418 _IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 419 flags=re.MULTILINE) 420 def _have_Iqxy(sources): 421 # type: (List[str]) -> bool 422 """ 423 Return true if any file defines Iqxy. 424 425 Note this is not a C parser, and so can be easily confused by 426 non-standard syntax. Also, it will incorrectly identify the following 427 as having Iqxy:: 428 429 /* 430 double Iqxy(qx, qy, ...) { ... fill this in later ... } 431 */ 432 433 If you want to comment out an Iqxy function, use // on the front of the 434 line instead. 435 """ 436 for code in sources: 437 if _IQXY_PATTERN.search(code): 438 return True 439 else: 440 return False 441 437 442 def make_source(model_info): 443 # type: (ModelInfo) -> str 438 444 """ 439 445 Generate the OpenCL/ctypes kernel from the module info. 440 446 441 Uses source files found in the given search path. 442 """ 443 if callable(model_info['Iq']): 444 return None 447 Uses source files found in the given search path. Returns None if this 448 is a pure python model, with no C source components. 449 """ 450 if callable(model_info.Iq): 451 raise ValueError("can't compile python model") 445 452 446 453 # TODO: need something other than volume to indicate dispersion parameters … … 453 460 # for computing volume even if we allow non-disperse volume parameters. 454 461 455 # Load template 456 global C_KERNEL_TEMPLATE 457 if C_KERNEL_TEMPLATE is None: 458 with open(C_KERNEL_TEMPLATE_PATH) as fid: 459 C_KERNEL_TEMPLATE = fid.read() 460 461 # Load additional sources 462 source = [open(f).read() for f in model_sources(model_info)] 463 464 # Prepare defines 465 defines = [] 466 partype = model_info['partype'] 467 pd_1d = partype['pd-1d'] 468 pd_2d = partype['pd-2d'] 469 fixed_1d = partype['fixed-1d'] 470 fixed_2d = partype['fixed-1d'] 471 472 iq_parameters = [p.name 473 for p in model_info['parameters'][2:] # skip scale, background 474 if p.name in set(fixed_1d + pd_1d)] 475 iqxy_parameters = [p.name 476 for p in model_info['parameters'][2:] # skip scale, background 477 if p.name in set(fixed_2d + pd_2d)] 478 volume_parameters = [p.name 479 for p in model_info['parameters'] 480 if p.type == 'volume'] 481 482 # Fill in defintions for volume parameters 483 if volume_parameters: 484 defines.append(('VOLUME_PARAMETERS', 485 ','.join(volume_parameters))) 486 defines.append(('VOLUME_WEIGHT_PRODUCT', 487 '*'.join(p + '_w' for p in volume_parameters))) 488 489 # Generate form_volume function from body only 490 if model_info['form_volume'] is not None: 491 if volume_parameters: 492 vol_par_decl = ', '.join('double ' + p for p in volume_parameters) 493 else: 494 vol_par_decl = 'void' 495 defines.append(('VOLUME_PARAMETER_DECLARATIONS', 496 vol_par_decl)) 497 fn = """\ 498 double form_volume(VOLUME_PARAMETER_DECLARATIONS); 499 double form_volume(VOLUME_PARAMETER_DECLARATIONS) { 500 %(body)s 501 } 502 """ % {'body':model_info['form_volume']} 503 source.append(fn) 504 505 # Fill in definitions for Iq parameters 506 defines.append(('IQ_KERNEL_NAME', model_info['name'] + '_Iq')) 507 defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters))) 508 if fixed_1d: 509 defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS', 510 ', \\\n '.join('const double %s' % p for p in fixed_1d))) 511 if pd_1d: 512 defines.append(('IQ_WEIGHT_PRODUCT', 513 '*'.join(p + '_w' for p in pd_1d))) 514 defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS', 515 ', \\\n '.join('const int N%s' % p for p in pd_1d))) 516 defines.append(('IQ_DISPERSION_LENGTH_SUM', 517 '+'.join('N' + p for p in pd_1d))) 518 open_loops, close_loops = build_polydispersity_loops(pd_1d) 519 defines.append(('IQ_OPEN_LOOPS', 520 open_loops.replace('\n', ' \\\n'))) 521 defines.append(('IQ_CLOSE_LOOPS', 522 close_loops.replace('\n', ' \\\n'))) 523 if model_info['Iq'] is not None: 524 defines.append(('IQ_PARAMETER_DECLARATIONS', 525 ', '.join('double ' + p for p in iq_parameters))) 526 fn = """\ 527 double Iq(double q, IQ_PARAMETER_DECLARATIONS); 528 double Iq(double q, IQ_PARAMETER_DECLARATIONS) { 529 %(body)s 530 } 531 """ % {'body':model_info['Iq']} 532 source.append(fn) 533 534 # Fill in definitions for Iqxy parameters 535 defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy')) 536 defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters))) 537 if fixed_2d: 538 defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS', 539 ', \\\n '.join('const double %s' % p for p in fixed_2d))) 540 if pd_2d: 541 defines.append(('IQXY_WEIGHT_PRODUCT', 542 '*'.join(p + '_w' for p in pd_2d))) 543 defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS', 544 ', \\\n '.join('const int N%s' % p for p in pd_2d))) 545 defines.append(('IQXY_DISPERSION_LENGTH_SUM', 546 '+'.join('N' + p for p in pd_2d))) 547 open_loops, close_loops = build_polydispersity_loops(pd_2d) 548 defines.append(('IQXY_OPEN_LOOPS', 549 open_loops.replace('\n', ' \\\n'))) 550 defines.append(('IQXY_CLOSE_LOOPS', 551 close_loops.replace('\n', ' \\\n'))) 552 if model_info['Iqxy'] is not None: 553 defines.append(('IQXY_PARAMETER_DECLARATIONS', 554 ', '.join('double ' + p for p in iqxy_parameters))) 555 fn = """\ 556 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS); 557 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) { 558 %(body)s 559 } 560 """ % {'body':model_info['Iqxy']} 561 source.append(fn) 562 563 # Need to know if we have a theta parameter for Iqxy; it is not there 564 # for the magnetic sphere model, for example, which has a magnetic 565 # orientation but no shape orientation. 566 if 'theta' in pd_2d: 567 defines.append(('IQXY_HAS_THETA', '1')) 568 569 #for d in defines: print(d) 570 defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 571 sources = '\n\n'.join(source) 572 return C_KERNEL_TEMPLATE % { 573 'DEFINES': defines, 574 'SOURCES': sources, 575 } 576 577 def categorize_parameters(pars): 578 """ 579 Build parameter categories out of the the parameter definitions. 580 581 Returns a dictionary of categories. 582 583 Note: these categories are subject to change, depending on the needs of 584 the UI and the needs of the kernel calling function. 585 586 The categories are as follows: 587 588 * *volume* list of volume parameter names 589 * *orientation* list of orientation parameters 590 * *magnetic* list of magnetic parameters 591 * *<empty string>* list of parameters that have no type info 592 593 Each parameter is in one and only one category. 594 595 The following derived categories are created: 596 597 * *fixed-1d* list of non-polydisperse parameters for 1D models 598 * *pd-1d* list of polydisperse parameters for 1D models 599 * *fixed-2d* list of non-polydisperse parameters for 2D models 600 * *pd-d2* list of polydisperse parameters for 2D models 601 """ 602 partype = { 603 'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], '': [], 604 'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [], 605 'pd-rel': set(), 606 } 607 608 for p in pars: 609 if p.type == 'volume': 610 partype['pd-1d'].append(p.name) 611 partype['pd-2d'].append(p.name) 612 partype['pd-rel'].add(p.name) 613 elif p.type == 'magnetic': 614 partype['fixed-2d'].append(p.name) 615 elif p.type == 'orientation': 616 partype['pd-2d'].append(p.name) 617 elif p.type in ('', 'sld'): 618 partype['fixed-1d'].append(p.name) 619 partype['fixed-2d'].append(p.name) 620 else: 621 raise ValueError("unknown parameter type %r" % p.type) 622 partype[p.type].append(p.name) 623 624 return partype 625 626 def process_parameters(model_info): 627 """ 628 Process parameter block, precalculating parameter details. 629 """ 630 # convert parameters into named tuples 631 for p in model_info['parameters']: 632 if p[4] == '' and (p[0].startswith('sld') or p[0].endswith('sld')): 633 p[4] = 'sld' 634 # TODO: make sure all models explicitly label their sld parameters 635 #raise ValueError("%s.%s needs to be explicitly set to type 'sld'" %(model_info['id'], p[0])) 636 637 pars = [Parameter(*p) for p in model_info['parameters']] 638 # Fill in the derived attributes 639 model_info['parameters'] = pars 640 partype = categorize_parameters(pars) 641 model_info['limits'] = dict((p.name, p.limits) for p in pars) 642 model_info['partype'] = partype 643 model_info['defaults'] = dict((p.name, p.default) for p in pars) 644 if model_info.get('demo', None) is None: 645 model_info['demo'] = model_info['defaults'] 646 model_info['has_2d'] = partype['orientation'] or partype['magnetic'] 647 462 partable = model_info.parameters 463 464 # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume. 465 # Note that scale and volume are not possible types. 466 467 # Load templates and user code 468 kernel_header = load_template('kernel_header.c') 469 kernel_code = load_template('kernel_iq.c') 470 user_code = [open(f).read() for f in model_sources(model_info)] 471 472 # Build initial sources 473 source = [kernel_header] + user_code 474 475 # Make parameters for q, qx, qy so that we can use them in declarations 476 q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 477 # Generate form_volume function, etc. from body only 478 if isinstance(model_info.form_volume, str): 479 pars = partable.form_volume_parameters 480 source.append(_gen_fn('form_volume', pars, model_info.form_volume)) 481 if isinstance(model_info.Iq, str): 482 pars = [q] + partable.iq_parameters 483 source.append(_gen_fn('Iq', pars, model_info.Iq)) 484 if isinstance(model_info.Iqxy, str): 485 pars = [qx, qy] + partable.iqxy_parameters 486 source.append(_gen_fn('Iqxy', pars, model_info.Iqxy)) 487 488 # Define the parameter table 489 source.append("#define PARAMETER_TABLE \\") 490 source.append("\\\n".join(p.as_definition() 491 for p in partable.kernel_parameters)) 492 493 # Define the function calls 494 if partable.form_volume_parameters: 495 refs = _call_pars("_v.", partable.form_volume_parameters) 496 call_volume = "#define CALL_VOLUME(_v) form_volume(%s)" % (",".join(refs)) 497 else: 498 # Model doesn't have volume. We could make the kernel run a little 499 # faster by not using/transferring the volume normalizations, but 500 # the ifdef's reduce readability more than is worthwhile. 501 call_volume = "#define CALL_VOLUME(v) 1.0" 502 source.append(call_volume) 503 504 refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 505 call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 506 if _have_Iqxy(user_code): 507 # Call 2D model 508 refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 509 call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 510 else: 511 # Call 1D model with sqrt(qx^2 + qy^2) 512 warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 513 # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 514 pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:] 515 call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 516 517 # Fill in definitions for numbers of parameters 518 source.append("#define MAX_PD %s"%partable.max_pd) 519 source.append("#define NPARS %d"%partable.npars) 520 521 # TODO: allow mixed python/opencl kernels? 522 523 # define the Iq kernel 524 source.append("#define KERNEL_NAME %s_Iq"%model_info.name) 525 source.append(call_iq) 526 source.append(kernel_code) 527 source.append("#undef CALL_IQ") 528 source.append("#undef KERNEL_NAME") 529 530 # define the Iqxy kernel from the same source with different #defines 531 source.append("#define KERNEL_NAME %s_Iqxy"%model_info.name) 532 source.append(call_iqxy) 533 source.append(kernel_code) 534 source.append("#undef CALL_IQ") 535 source.append("#undef KERNEL_NAME") 536 537 return '\n'.join(source) 648 538 649 539 def load_kernel_module(model_name): 540 # type: (str) -> module 650 541 if model_name.endswith('.py'): 651 542 kernel_module = load_custom_kernel_module(model_name) … … 657 548 658 549 659 def make_model_info(kernel_module):660 """661 Interpret the model definition file, categorizing the parameters.662 663 The module can be loaded with a normal python import statement if you664 know which module you need, or with __import__('sasmodels.model.'+name)665 if the name is in a string.666 667 The *model_info* structure contains the following fields:668 669 * *id* is the id of the kernel670 * *name* is the display name of the kernel671 * *filename* is the full path to the module defining the file (if any)672 * *title* is a short description of the kernel673 * *description* is a long description of the kernel (this doesn't seem674 very useful since the Help button on the model page brings you directly675 to the documentation page)676 * *docs* is the docstring from the module. Use :func:`make_doc` to677 * *category* specifies the model location in the docs678 * *parameters* is the model parameter table679 * *single* is True if the model allows single precision680 * *structure_factor* is True if the model is useable in a product681 * *defaults* is the *{parameter: value}* table built from the parameter682 description table.683 * *limits* is the *{parameter: [min, max]}* table built from the684 parameter description table.685 * *partypes* categorizes the model parameters. See686 :func:`categorize_parameters` for details.687 * *demo* contains the *{parameter: value}* map used in compare (and maybe688 for the demo plot, if plots aren't set up to use the default values).689 If *demo* is not given in the file, then the default values will be used.690 * *tests* is a set of tests that must pass691 * *source* is the list of library files to include in the C model build692 * *Iq*, *Iqxy*, *form_volume*, *ER*, *VR* and *sesans* are python functions693 implementing the kernel for the module, or None if they are not694 defined in python695 * *composition* is None if the model is independent, otherwise it is a696 tuple with composition type ('product' or 'mixture') and a list of697 *model_info* blocks for the composition objects. This allows us to698 build complete product and mixture models from just the info.699 700 """701 parameters = COMMON_PARAMETERS + kernel_module.parameters702 filename = abspath(kernel_module.__file__)703 kernel_id = splitext(basename(filename))[0]704 name = getattr(kernel_module, 'name', None)705 if name is None:706 name = " ".join(w.capitalize() for w in kernel_id.split('_'))707 model_info = dict(708 id=kernel_id, # string used to load the kernel709 filename=abspath(kernel_module.__file__),710 name=name,711 title=kernel_module.title,712 description=kernel_module.description,713 parameters=parameters,714 composition=None,715 docs=kernel_module.__doc__,716 category=getattr(kernel_module, 'category', None),717 single=getattr(kernel_module, 'single', True),718 structure_factor=getattr(kernel_module, 'structure_factor', False),719 control=getattr(kernel_module, 'control', None),720 demo=getattr(kernel_module, 'demo', None),721 source=getattr(kernel_module, 'source', []),722 tests=getattr(kernel_module, 'tests', []),723 )724 process_parameters(model_info)725 # Check for optional functions726 functions = "ER VR form_volume Iq Iqxy shape sesans".split()727 model_info.update((k, getattr(kernel_module, k, None)) for k in functions)728 return model_info729 550 730 551 section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z' 731 552 %re.escape(string.punctuation)) 732 553 def _convert_section_titles_to_boldface(lines): 554 # type: (Sequence[str]) -> Iterator[str] 733 555 """ 734 556 Do the actual work of identifying and converting section headings. … … 752 574 753 575 def convert_section_titles_to_boldface(s): 576 # type: (str) -> str 754 577 """ 755 578 Use explicit bold-face rather than section headings so that the table of … … 762 585 763 586 def make_doc(model_info): 587 # type: (ModelInfo) -> str 764 588 """ 765 589 Return the documentation for the model. … … 767 591 Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale." 768 592 Sq_units = "The returned value is a dimensionless structure factor, $S(q)$." 769 docs = convert_section_titles_to_boldface(model_info['docs']) 770 subst = dict(id=model_info['id'].replace('_', '-'), 771 name=model_info['name'], 772 title=model_info['title'], 773 parameters=make_partable(model_info['parameters']), 774 returns=Sq_units if model_info['structure_factor'] else Iq_units, 593 docs = convert_section_titles_to_boldface(model_info.docs) 594 pars = make_partable(model_info.parameters.COMMON 595 + model_info.parameters.kernel_parameters) 596 subst = dict(id=model_info.id.replace('_', '-'), 597 name=model_info.name, 598 title=model_info.title, 599 parameters=pars, 600 returns=Sq_units if model_info.structure_factor else Iq_units, 775 601 docs=docs) 776 602 return DOC_HEADER % subst … … 778 604 779 605 def demo_time(): 606 # type: () -> None 780 607 """ 781 608 Show how long it takes to process a model. 782 609 """ 610 import datetime 611 from .modelinfo import make_model_info 783 612 from .models import cylinder 784 import datetime 613 785 614 tic = datetime.datetime.now() 786 615 make_source(make_model_info(cylinder)) … … 789 618 790 619 def main(): 620 # type: () -> None 791 621 """ 792 622 Program which prints the source produced by the model. 793 623 """ 624 import sys 625 from .modelinfo import make_model_info 626 794 627 if len(sys.argv) <= 1: 795 628 print("usage: python -m sasmodels.generate modelname") -
sasmodels/kernelcl.py
r4d76711 ra5b8477 48 48 harmless, albeit annoying. 49 49 """ 50 from __future__ import print_function 51 50 52 import os 51 53 import warnings 52 54 53 import numpy as np 55 import numpy as np # type: ignore 54 56 55 57 try: 56 import pyopencl as cl 58 raise NotImplementedError("OpenCL not yet implemented for new kernel template") 59 import pyopencl as cl # type: ignore 57 60 # Ask OpenCL for the default context so that we know that one exists 58 61 cl.create_some_context(interactive=False) … … 61 64 raise RuntimeError("OpenCL not available") 62 65 63 from pyopencl import mem_flags as mf 64 from pyopencl.characterize import get_fast_inaccurate_build_options 66 from pyopencl import mem_flags as mf # type: ignore 67 from pyopencl.characterize import get_fast_inaccurate_build_options # type: ignore 65 68 66 69 from . import generate 70 from .kernel import KernelModel, Kernel 71 72 try: 73 from typing import Tuple, Callable, Any 74 from .modelinfo import ModelInfo 75 from .details import CallDetails 76 except ImportError: 77 pass 67 78 68 79 # The max loops number is limited by the amount of local memory available … … 73 84 # of polydisperse parameters. 74 85 MAX_LOOPS = 2048 86 87 88 # Pragmas for enable OpenCL features. Be sure to protect them so that they 89 # still compile even if OpenCL is not present. 90 _F16_PRAGMA = """\ 91 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 92 # pragma OPENCL EXTENSION cl_khr_fp16: enable 93 #endif 94 """ 95 96 _F64_PRAGMA = """\ 97 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 98 # pragma OPENCL EXTENSION cl_khr_fp64: enable 99 #endif 100 """ 75 101 76 102 … … 142 168 raise RuntimeError("%s not supported for devices"%dtype) 143 169 144 source = generate.convert_type(source, dtype) 170 source_list = [generate.convert_type(source, dtype)] 171 172 if dtype == generate.F16: 173 source_list.insert(0, _F16_PRAGMA) 174 elif dtype == generate.F64: 175 source_list.insert(0, _F64_PRAGMA) 176 145 177 # Note: USE_SINCOS makes the intel cpu slower under opencl 146 178 if context.devices[0].type == cl.device_type.GPU: 147 source = "#define USE_SINCOS\n" + source179 source_list.insert(0, "#define USE_SINCOS\n") 148 180 options = (get_fast_inaccurate_build_options(context.devices[0]) 149 181 if fast else []) 182 source = "\n".join(source_list) 150 183 program = cl.Program(context, source).build(options=options) 151 184 return program … … 220 253 key = "%s-%s-%s"%(name, dtype, fast) 221 254 if key not in self.compiled: 222 #print("compiling",name)255 print("compiling",name) 223 256 dtype = np.dtype(dtype) 224 program = compile_model(self.get_context(dtype), source, dtype, fast) 257 program = compile_model(self.get_context(dtype), 258 str(source), dtype, fast) 225 259 self.compiled[key] = program 226 260 return self.compiled[key] … … 285 319 286 320 287 class GpuModel( object):321 class GpuModel(KernelModel): 288 322 """ 289 323 GPU wrapper for a single model. … … 317 351 if self.program is None: 318 352 compiler = environment().compile_program 319 self.program = compiler(self.info ['name'], self.source, self.dtype,320 self. fast)353 self.program = compiler(self.info.name, self.source, 354 self.dtype, self.fast) 321 355 is_2d = len(q_vectors) == 2 322 356 kernel_name = generate.kernel_name(self.info, is_2d) … … 329 363 """ 330 364 if self.program is not None: 331 environment().release_program(self.info ['name'])365 environment().release_program(self.info.name) 332 366 self.program = None 333 367 … … 365 399 # at this point, so instead using 32, which is good on the set of 366 400 # architectures tested so far. 367 self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 401 if self.is_2d: 402 # Note: 17 rather than 15 because results is 2 elements 403 # longer than input. 404 width = ((self.nq+17)//16)*16 405 self.q = np.empty((width, 2), dtype=dtype) 406 self.q[:self.nq, 0] = q_vectors[0] 407 self.q[:self.nq, 1] = q_vectors[1] 408 else: 409 # Note: 33 rather than 31 because results is 2 elements 410 # longer than input. 411 width = ((self.nq+33)//32)*32 412 self.q = np.empty(width, dtype=dtype) 413 self.q[:self.nq] = q_vectors[0] 414 self.global_size = [self.q.shape[0]] 368 415 context = env.get_context(self.dtype) 369 self.global_size = [self.q_vectors[0].size]370 416 #print("creating inputs of size", self.global_size) 371 self.q_buffers = [ 372 cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 373 for q in self.q_vectors 374 ] 417 self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 418 hostbuf=self.q) 375 419 376 420 def release(self): … … 378 422 Free the memory. 379 423 """ 380 for b in self.q_buffers:381 b.release()382 self.q_buffers = []424 if self.q is not None: 425 self.q.release() 426 self.q = None 383 427 384 428 def __del__(self): 385 429 self.release() 386 430 387 class GpuKernel( object):431 class GpuKernel(Kernel): 388 432 """ 389 433 Callable SAS kernel. … … 405 449 Call :meth:`release` when done with the kernel instance. 406 450 """ 407 def __init__(self, kernel, model_info, q_vectors, dtype): 408 q_input = GpuInput(q_vectors, dtype) 451 def __init__(self, kernel, model_info, q_vectors): 452 # type: (KernelModel, ModelInfo, List[np.ndarray]) -> None 453 max_pd = model_info.parameters.max_pd 454 npars = len(model_info.parameters.kernel_parameters)-2 455 q_input = GpuInput(q_vectors, kernel.dtype) 409 456 self.kernel = kernel 410 457 self.info = model_info 411 self.res = np.empty(q_input.nq, q_input.dtype) 412 dim = '2d' if q_input.is_2d else '1d' 413 self.fixed_pars = model_info['partype']['fixed-' + dim] 414 self.pd_pars = model_info['partype']['pd-' + dim] 458 self.dtype = kernel.dtype 459 self.dim = '2d' if q_input.is_2d else '1d' 460 self.pd_stop_index = 4*max_pd-1 461 # plus three for the normalization values 462 self.result = np.empty(q_input.nq+3, q_input.dtype) 415 463 416 464 # Inputs and outputs for each kernel call 417 465 # Note: res may be shorter than res_b if global_size != nq 418 466 env = environment() 419 self.queue = env.get_queue(dtype) 420 self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 421 2 * MAX_LOOPS * q_input.dtype.itemsize) 422 self.res_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 423 q_input.global_size[0] * q_input.dtype.itemsize) 424 self.q_input = q_input 425 426 self._need_release = [self.loops_b, self.res_b, self.q_input] 427 428 def __call__(self, fixed_pars, pd_pars, cutoff): 467 self.queue = env.get_queue(kernel.dtype) 468 469 # details is int32 data, padded to an 8 integer boundary 470 size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8 471 self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 472 q_input.global_size[0] * kernel.dtype.itemsize) 473 self.q_input = q_input # allocated by GpuInput above 474 475 self._need_release = [ self.result_b, self.q_input ] 476 477 def __call__(self, call_details, weights, values, cutoff): 478 # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 429 479 real = (np.float32 if self.q_input.dtype == generate.F32 430 480 else np.float64 if self.q_input.dtype == generate.F64 431 481 else np.float16 if self.q_input.dtype == generate.F16 432 482 else np.float32) # will never get here, so use np.float32 433 434 #print "pars", fixed_pars, pd_pars 435 res_bi = self.res_b 436 nq = np.uint32(self.q_input.nq) 437 if pd_pars: 438 cutoff = real(cutoff) 439 loops_N = [np.uint32(len(p[0])) for p in pd_pars] 440 loops = np.hstack(pd_pars) \ 441 if pd_pars else np.empty(0, dtype=self.q_input.dtype) 442 loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 443 #print("loops",Nloops, loops) 444 445 #import sys; print("opencl eval",pars) 446 #print("opencl eval",pars) 447 if len(loops) > 2 * MAX_LOOPS: 448 raise ValueError("too many polydispersity points") 449 450 loops_bi = self.loops_b 451 cl.enqueue_copy(self.queue, loops_bi, loops) 452 loops_l = cl.LocalMemory(len(loops.data)) 453 #ctx = environment().context 454 #loops_bi = cl.Buffer(ctx, mf.READ_ONLY|mf.COPY_HOST_PTR, hostbuf=loops) 455 dispersed = [loops_bi, loops_l, cutoff] + loops_N 456 else: 457 dispersed = [] 458 fixed = [real(p) for p in fixed_pars] 459 args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed 483 assert call_details.dtype == np.int32 484 assert weights.dtype == real and values.dtype == real 485 486 context = self.queue.context 487 details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 488 hostbuf=call_details) 489 weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 490 hostbuf=weights) 491 values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 492 hostbuf=values) 493 494 start, stop = 0, self.details[self.pd_stop_index] 495 args = [ 496 np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop), 497 self.details_b, self.weights_b, self.values_b, 498 self.q_input.q_b, self.result_b, real(cutoff), 499 ] 460 500 self.kernel(self.queue, self.q_input.global_size, None, *args) 461 cl.enqueue_copy(self.queue, self.res, res_bi) 462 463 return self.res 501 cl.enqueue_copy(self.queue, self.result, self.result_b) 502 [v.release() for v in (details_b, weights_b, values_b)] 503 504 return self.result[:self.nq] 464 505 465 506 def release(self): -
sasmodels/kerneldll.py
r4d76711 ra5b8477 49 49 import os 50 50 import tempfile 51 import ctypes as ct 52 from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float 53 import _ctypes 54 55 import numpy as np 51 import ctypes as ct # type: ignore 52 from ctypes import c_void_p, c_int32, c_longdouble, c_double, c_float # type: ignore 53 54 import numpy as np # type: ignore 56 55 57 56 from . import generate 58 from .kernelpy import PyInput, PyModel 57 from .kernel import KernelModel, Kernel 58 from .kernelpy import PyInput 59 59 from .exception import annotate_exception 60 from .generate import F16, F32, F64 61 62 try: 63 from typing import Tuple, Callable, Any 64 from .modelinfo import ModelInfo 65 from .details import CallDetails 66 except ImportError: 67 pass 60 68 61 69 # Compiler platform details … … 81 89 COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 82 90 if "SAS_OPENMP" in os.environ: 83 COMPILE = COMPILE +" -fopenmp"91 COMPILE += " -fopenmp" 84 92 else: 85 93 COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" … … 90 98 91 99 92 def dll_path(model_info, dtype="double"): 93 """ 94 Path to the compiled model defined by *model_info*. 95 """ 96 from os.path import join as joinpath, split as splitpath, splitext 97 basename = splitext(splitpath(model_info['filename'])[1])[0] 98 if np.dtype(dtype) == generate.F32: 99 basename += "32" 100 elif np.dtype(dtype) == generate.F64: 101 basename += "64" 102 else: 103 basename += "128" 104 return joinpath(DLL_PATH, basename+'.so') 105 106 107 def make_dll(source, model_info, dtype="double"): 108 """ 109 Load the compiled model defined by *kernel_module*. 110 111 Recompile if any files are newer than the model file. 100 def dll_name(model_info, dtype): 101 # type: (ModelInfo, np.dtype) -> str 102 """ 103 Name of the dll containing the model. This is the base file name without 104 any path or extension, with a form such as 'sas_sphere32'. 105 """ 106 bits = 8*dtype.itemsize 107 return "sas_%s%d"%(model_info.id, bits) 108 109 110 def dll_path(model_info, dtype): 111 # type: (ModelInfo, np.dtype) -> str 112 """ 113 Complete path to the dll for the model. Note that the dll may not 114 exist yet if it hasn't been compiled. 115 """ 116 return os.path.join(DLL_PATH, dll_name(model_info, dtype)+".so") 117 118 119 def make_dll(source, model_info, dtype=F64): 120 # type: (str, ModelInfo, np.dtype) -> str 121 """ 122 Returns the path to the compiled model defined by *kernel_module*. 123 124 If the model has not been compiled, or if the source file(s) are newer 125 than the dll, then *make_dll* will compile the model before returning. 126 This routine does not load the resulting dll. 112 127 113 128 *dtype* is a numpy floating point precision specifier indicating whether 114 the model should be single or double precision. The default is double115 precision.116 117 The DLL is not loaded until the kernel is called so models can118 be defined without using too many resources.129 the model should be single, double or long double precision. The default 130 is double precision, *np.dtype('d')*. 131 132 Set *sasmodels.ALLOW_SINGLE_PRECISION_DLLS* to False if single precision 133 models are not allowed as DLLs. 119 134 120 135 Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 121 136 The default is the system temporary directory. 122 123 Set *sasmodels.ALLOW_SINGLE_PRECISION_DLLS* to True if single precision 124 models are allowed as DLLs. 125 """ 126 if callable(model_info.get('Iq', None)): 127 return PyModel(model_info) 128 129 dtype = np.dtype(dtype) 130 if dtype == generate.F16: 137 """ 138 if dtype == F16: 131 139 raise ValueError("16 bit floats not supported") 132 if dtype == generate.F32 and not ALLOW_SINGLE_PRECISION_DLLS: 133 dtype = generate.F64 # Force 64-bit dll 134 135 if dtype == generate.F32: # 32-bit dll 136 tempfile_prefix = 'sas_' + model_info['name'] + '32_' 137 elif dtype == generate.F64: 138 tempfile_prefix = 'sas_' + model_info['name'] + '64_' 139 else: 140 tempfile_prefix = 'sas_' + model_info['name'] + '128_' 141 142 source = generate.convert_type(source, dtype) 143 source_files = generate.model_sources(model_info) + [model_info['filename']] 140 if dtype == F32 and not ALLOW_SINGLE_PRECISION_DLLS: 141 dtype = F64 # Force 64-bit dll 142 # Note: dtype may be F128 for long double precision 143 144 newest = generate.timestamp(model_info) 144 145 dll = dll_path(model_info, dtype) 145 newest = max(os.path.getmtime(f) for f in source_files)146 146 if not os.path.exists(dll) or os.path.getmtime(dll) < newest: 147 # Replace with a proper temp file 148 fid, filename = tempfile.mkstemp(suffix=".c", prefix=tempfile_prefix) 147 basename = dll_name(model_info, dtype) + "_" 148 fid, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 149 source = generate.convert_type(source, dtype) 149 150 os.fdopen(fid, "w").write(source) 150 151 command = COMPILE%{"source":filename, "output":dll} … … 160 161 161 162 162 def load_dll(source, model_info, dtype="double"): 163 def load_dll(source, model_info, dtype=F64): 164 # type: (str, ModelInfo, np.dtype) -> "DllModel" 163 165 """ 164 166 Create and load a dll corresponding to the source, info pair returned … … 172 174 173 175 174 IQ_ARGS = [c_void_p, c_void_p, c_int] 175 IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int] 176 177 class DllModel(object): 176 class DllModel(KernelModel): 178 177 """ 179 178 ctypes wrapper for a single model. … … 191 190 192 191 def __init__(self, dllpath, model_info, dtype=generate.F32): 192 # type: (str, ModelInfo, np.dtype) -> None 193 193 self.info = model_info 194 194 self.dllpath = dllpath 195 self. dll = None195 self._dll = None # type: ct.CDLL 196 196 self.dtype = np.dtype(dtype) 197 197 198 198 def _load_dll(self): 199 Nfixed1d = len(self.info['partype']['fixed-1d']) 200 Nfixed2d = len(self.info['partype']['fixed-2d']) 201 Npd1d = len(self.info['partype']['pd-1d']) 202 Npd2d = len(self.info['partype']['pd-2d']) 203 199 # type: () -> None 204 200 #print("dll", self.dllpath) 205 201 try: 206 self. dll = ct.CDLL(self.dllpath)202 self._dll = ct.CDLL(self.dllpath) 207 203 except: 208 204 annotate_exception("while loading "+self.dllpath) … … 212 208 else c_double if self.dtype == generate.F64 213 209 else c_longdouble) 214 pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else [] 215 pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else [] 216 self.Iq = self.dll[generate.kernel_name(self.info, False)] 217 self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d 218 219 self.Iqxy = self.dll[generate.kernel_name(self.info, True)] 220 self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [fp]*Nfixed2d 221 222 self.release() 210 211 # int, int, int, int*, double*, double*, double*, double*, double*, double 212 argtypes = [c_int32]*3 + [c_void_p]*5 + [fp] 213 self._Iq = self._dll[generate.kernel_name(self.info, is_2d=False)] 214 self._Iqxy = self._dll[generate.kernel_name(self.info, is_2d=True)] 215 self._Iq.argtypes = argtypes 216 self._Iqxy.argtypes = argtypes 223 217 224 218 def __getstate__(self): 219 # type: () -> Tuple[ModelInfo, str] 225 220 return self.info, self.dllpath 226 221 227 222 def __setstate__(self, state): 223 # type: (Tuple[ModelInfo, str]) -> None 228 224 self.info, self.dllpath = state 229 self. dll = None225 self._dll = None 230 226 231 227 def make_kernel(self, q_vectors): 228 # type: (List[np.ndarray]) -> DllKernel 232 229 q_input = PyInput(q_vectors, self.dtype) 233 if self.dll is None: self._load_dll() 234 kernel = self.Iqxy if q_input.is_2d else self.Iq 230 # Note: pickle not supported for DllKernel 231 if self._dll is None: 232 self._load_dll() 233 kernel = self._Iqxy if q_input.is_2d else self._Iq 235 234 return DllKernel(kernel, self.info, q_input) 236 235 237 236 def release(self): 237 # type: () -> None 238 238 """ 239 239 Release any resources associated with the model. … … 244 244 libHandle = dll._handle 245 245 #libHandle = ct.c_void_p(dll._handle) 246 del dll, self. dll247 self. dll = None246 del dll, self._dll 247 self._dll = None 248 248 #_ctypes.FreeLibrary(libHandle) 249 249 ct.windll.kernel32.FreeLibrary(libHandle) … … 252 252 253 253 254 class DllKernel( object):254 class DllKernel(Kernel): 255 255 """ 256 256 Callable SAS kernel. … … 272 272 """ 273 273 def __init__(self, kernel, model_info, q_input): 274 # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 275 self.kernel = kernel 274 276 self.info = model_info 275 277 self.q_input = q_input 276 self.kernel = kernel 277 self.res = np.empty(q_input.nq, q_input.dtype) 278 dim = '2d' if q_input.is_2d else '1d' 279 self.fixed_pars = model_info['partype']['fixed-' + dim] 280 self.pd_pars = model_info['partype']['pd-' + dim] 281 282 # In dll kernel, but not in opencl kernel 283 self.p_res = self.res.ctypes.data 284 285 def __call__(self, fixed_pars, pd_pars, cutoff): 286 real = (np.float32 if self.q_input.dtype == generate.F32 287 else np.float64 if self.q_input.dtype == generate.F64 288 else np.float128) 289 290 nq = c_int(self.q_input.nq) 291 if pd_pars: 292 cutoff = real(cutoff) 293 loops_N = [np.uint32(len(p[0])) for p in pd_pars] 294 loops = np.hstack(pd_pars) 295 loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 296 p_loops = loops.ctypes.data 297 dispersed = [p_loops, cutoff] + loops_N 298 else: 299 dispersed = [] 300 fixed = [real(p) for p in fixed_pars] 301 args = self.q_input.q_pointers + [self.p_res, nq] + dispersed + fixed 302 #print(pars) 303 self.kernel(*args) 304 305 return self.res 278 self.dtype = q_input.dtype 279 self.dim = '2d' if q_input.is_2d else '1d' 280 self.result = np.empty(q_input.nq+1, q_input.dtype) 281 self.real = (np.float32 if self.q_input.dtype == generate.F32 282 else np.float64 if self.q_input.dtype == generate.F64 283 else np.float128) 284 285 def __call__(self, call_details, weights, values, cutoff): 286 # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 287 288 #print("in kerneldll") 289 #print("weights", weights) 290 #print("values", values) 291 start, stop = 0, call_details.total_pd 292 args = [ 293 self.q_input.nq, # nq 294 start, # pd_start 295 stop, # pd_stop pd_stride[MAX_PD] 296 call_details.ctypes.data, # problem 297 weights.ctypes.data, # weights 298 values.ctypes.data, #pars 299 self.q_input.q.ctypes.data, #q 300 self.result.ctypes.data, # results 301 self.real(cutoff), # cutoff 302 ] 303 self.kernel(*args) # type: ignore 304 return self.result[:-1] 306 305 307 306 def release(self): 307 # type: () -> None 308 308 """ 309 309 Release any resources associated with the kernel. 310 310 """ 311 pass311 self.q_input.release() -
sasmodels/kernelpy.py
r4d76711 r7ae2b7f 7 7 :class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`. 8 8 """ 9 import numpy as np 10 from numpy import pi, cos 11 9 import numpy as np # type: ignore 10 from numpy import pi, cos #type: ignore 11 12 from . import details 12 13 from .generate import F64 13 14 class PyModel(object): 14 from .kernel import KernelModel, Kernel 15 16 try: 17 from typing import Union, Callable 18 except: 19 pass 20 else: 21 DType = Union[None, str, np.dtype] 22 23 class PyModel(KernelModel): 15 24 """ 16 25 Wrapper for pure python models. 17 26 """ 18 27 def __init__(self, model_info): 28 # Make sure Iq and Iqxy are available and vectorized 29 _create_default_functions(model_info) 19 30 self.info = model_info 20 31 21 32 def make_kernel(self, q_vectors): 22 33 q_input = PyInput(q_vectors, dtype=F64) 23 kernel = self.info ['Iqxy'] if q_input.is_2d else self.info['Iq']34 kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 24 35 return PyKernel(kernel, self.info, q_input) 25 36 … … 53 64 self.dtype = dtype 54 65 self.is_2d = (len(q_vectors) == 2) 55 self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors] 56 self.q_pointers = [q.ctypes.data for q in self.q_vectors] 66 if self.is_2d: 67 self.q = np.empty((self.nq, 2), dtype=dtype) 68 self.q[:, 0] = q_vectors[0] 69 self.q[:, 1] = q_vectors[1] 70 else: 71 self.q = np.empty(self.nq, dtype=dtype) 72 self.q[:self.nq] = q_vectors[0] 57 73 58 74 def release(self): … … 60 76 Free resources associated with the model inputs. 61 77 """ 62 self.q _vectors = []63 64 class PyKernel( object):78 self.q = None 79 80 class PyKernel(Kernel): 65 81 """ 66 82 Callable SAS kernel. … … 82 98 """ 83 99 def __init__(self, kernel, model_info, q_input): 100 self.dtype = np.dtype('d') 84 101 self.info = model_info 85 102 self.q_input = q_input 86 103 self.res = np.empty(q_input.nq, q_input.dtype) 87 dim = '2d' if q_input.is_2d else '1d' 88 # Loop over q unless user promises that the kernel is vectorized by 89 # taggining it with vectorized=True 90 if not getattr(kernel, 'vectorized', False): 91 if dim == '2d': 92 def vector_kernel(qx, qy, *args): 93 """ 94 Vectorized 2D kernel. 95 """ 96 return np.array([kernel(qxi, qyi, *args) 97 for qxi, qyi in zip(qx, qy)]) 104 self.kernel = kernel 105 self.dim = '2d' if q_input.is_2d else '1d' 106 107 partable = model_info.parameters 108 kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 109 else partable.iq_parameters) 110 volume_parameters = partable.form_volume_parameters 111 112 # Create an array to hold the parameter values. There will be a 113 # single array whose values are updated as the calculator goes 114 # through the loop. Arguments to the kernel and volume functions 115 # will use views into this vector, relying on the fact that a 116 # an array of no dimensions acts like a scalar. 117 parameter_vector = np.empty(len(partable.call_parameters)-2, 'd') 118 119 # Create views into the array to hold the arguments 120 offset = 0 121 kernel_args, volume_args = [], [] 122 for p in partable.kernel_parameters: 123 if p.length == 1: 124 # Scalar values are length 1 vectors with no dimensions. 125 v = parameter_vector[offset:offset+1].reshape(()) 98 126 else: 99 def vector_kernel(q, *args): 100 """ 101 Vectorized 1D kernel. 102 """ 103 return np.array([kernel(qi, *args) 104 for qi in q]) 105 self.kernel = vector_kernel 127 # Vector values are simple views. 128 v = parameter_vector[offset:offset+p.length] 129 offset += p.length 130 if p in kernel_parameters: 131 kernel_args.append(v) 132 if p in volume_parameters: 133 volume_args.append(v) 134 135 # Hold on to the parameter vector so we can use it to call kernel later. 136 # This may also be required to preserve the views into the vector. 137 self._parameter_vector = parameter_vector 138 139 # Generate a closure which calls the kernel with the views into the 140 # parameter array. 141 if q_input.is_2d: 142 form = model_info.Iqxy 143 qx, qy = q_input.q[:,0], q_input.q[:,1] 144 self._form = lambda: form(qx, qy, *kernel_args) 106 145 else: 107 self.kernel = kernel 108 fixed_pars = model_info['partype']['fixed-' + dim] 109 pd_pars = model_info['partype']['pd-' + dim] 110 vol_pars = model_info['partype']['volume'] 111 112 # First two fixed pars are scale and background 113 pars = [p.name for p in model_info['parameters'][2:]] 114 offset = len(self.q_input.q_vectors) 115 self.args = self.q_input.q_vectors + [None] * len(pars) 116 self.fixed_index = np.array([pars.index(p) + offset 117 for p in fixed_pars[2:]]) 118 self.pd_index = np.array([pars.index(p) + offset 119 for p in pd_pars]) 120 self.vol_index = np.array([pars.index(p) + offset 121 for p in vol_pars]) 122 try: self.theta_index = pars.index('theta') + offset 123 except ValueError: self.theta_index = -1 124 125 # Caller needs fixed_pars and pd_pars 126 self.fixed_pars = fixed_pars 127 self.pd_pars = pd_pars 128 129 def __call__(self, fixed, pd, cutoff=1e-5): 130 #print("fixed",fixed) 131 #print("pd", pd) 132 args = self.args[:] # grab a copy of the args 133 form, form_volume = self.kernel, self.info['form_volume'] 134 # First two fixed 135 scale, background = fixed[:2] 136 for index, value in zip(self.fixed_index, fixed[2:]): 137 args[index] = float(value) 138 res = _loops(form, form_volume, cutoff, scale, background, args, 139 pd, self.pd_index, self.vol_index, self.theta_index) 140 146 form = model_info.Iq 147 q = q_input.q 148 self._form = lambda: form(q, *kernel_args) 149 150 # Generate a closure which calls the form_volume if it exists. 151 form_volume = model_info.form_volume 152 self._volume = ((lambda: form_volume(*volume_args)) if form_volume 153 else (lambda: 1.0)) 154 155 def __call__(self, call_details, weights, values, cutoff): 156 assert isinstance(call_details, details.CallDetails) 157 res = _loops(self._parameter_vector, self._form, self._volume, 158 self.q_input.nq, call_details, weights, values, cutoff) 141 159 return res 142 160 … … 145 163 Free resources associated with the kernel. 146 164 """ 165 self.q_input.release() 147 166 self.q_input = None 148 167 149 def _loops(form, form_volume, cutoff, scale, background, 150 args, pd, pd_index, vol_index, theta_index): 151 """ 152 *form* is the name of the form function, which should be vectorized over 153 q, but otherwise have an interface like the opencl kernels, with the 154 q parameters first followed by the individual parameters in the order 155 given in model.parameters (see :mod:`sasmodels.generate`). 156 157 *form_volume* calculates the volume of the shape. *vol_index* gives 158 the list of volume parameters 159 160 *cutoff* ignores the corners of the dispersion hypercube 161 162 *scale*, *background* multiplies the resulting form and adds an offset 163 164 *args* is the prepopulated set of arguments to the form function, starting 165 with the q vectors, and including slots for all the parameters. The 166 values for the parameters will be substituted with values from the 167 dispersion functions. 168 169 *pd* is the list of dispersion parameters 170 171 *pd_index* are the indices of the dispersion parameters in *args* 172 173 *vol_index* are the indices of the volume parameters in *args* 174 175 *theta_index* is the index of the theta parameter for the sasview 176 spherical correction, or -1 if there is no angular dispersion 177 """ 178 168 def _loops(parameters, form, form_volume, nq, call_details, 169 weights, values, cutoff): 170 # type: (np.ndarray, Callable[[], np.ndarray], Callable[[], float], int, details.CallDetails, np.ndarray, np.ndarray, float) -> None 179 171 ################################################################ 180 172 # # … … 186 178 # # 187 179 ################################################################ 188 189 weight = np.empty(len(pd), 'd') 190 if weight.size > 0: 191 # weight vector, to be populated by polydispersity loops 192 193 # identify which pd parameters are volume parameters 194 vol_weight_index = np.array([(index in vol_index) for index in pd_index]) 195 196 # Sort parameters in decreasing order of pd length 197 Npd = np.array([len(pdi[0]) for pdi in pd], 'i') 198 order = np.argsort(Npd)[::-1] 199 stride = np.cumprod(Npd[order]) 200 pd = [pd[index] for index in order] 201 pd_index = pd_index[order] 202 vol_weight_index = vol_weight_index[order] 203 204 fast_value = pd[0][0] 205 fast_weight = pd[0][1] 180 parameters[:] = values[call_details.par_offset] 181 scale, background = values[0], values[1] 182 if call_details.num_active == 0: 183 norm = float(form_volume()) 184 if norm > 0.0: 185 return (scale/norm)*form() + background 186 else: 187 return np.ones(nq, 'd')*background 188 189 partial_weight = np.NaN 190 spherical_correction = 1.0 191 pd_stride = call_details.pd_stride[:call_details.num_active] 192 pd_length = call_details.pd_length[:call_details.num_active] 193 pd_offset = call_details.pd_offset[:call_details.num_active] 194 pd_index = np.empty_like(pd_offset) 195 offset = np.empty_like(call_details.par_offset) 196 theta = call_details.theta_par 197 fast_length = pd_length[0] 198 pd_index[0] = fast_length 199 total = np.zeros(nq, 'd') 200 norm = 0.0 201 for loop_index in range(call_details.total_pd): 202 # update polydispersity parameter values 203 if pd_index[0] == fast_length: 204 pd_index[:] = (loop_index/pd_stride)%pd_length 205 partial_weight = np.prod(weights[pd_offset+pd_index][1:]) 206 for k in range(call_details.num_coord): 207 par = call_details.par_coord[k] 208 coord = call_details.pd_coord[k] 209 this_offset = call_details.par_offset[par] 210 block_size = 1 211 for bit in range(len(pd_offset)): 212 if coord&1: 213 this_offset += block_size * pd_index[bit] 214 block_size *= pd_length[bit] 215 coord >>= 1 216 if coord == 0: break 217 offset[par] = this_offset 218 parameters[par] = values[this_offset] 219 if par == theta and not (call_details.par_coord[k]&1): 220 spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 221 for k in range(call_details.num_coord): 222 if call_details.pd_coord[k]&1: 223 #par = call_details.par_coord[k] 224 parameters[par] = values[offset[par]] 225 #print "par",par,offset[par],parameters[par+2] 226 offset[par] += 1 227 if par == theta: 228 spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 229 230 weight = partial_weight * weights[pd_offset[0] + pd_index[0]] 231 pd_index[0] += 1 232 if weight > cutoff: 233 # Call the scattering function 234 # Assume that NaNs are only generated if the parameters are bad; 235 # exclude all q for that NaN. Even better would be to have an 236 # INVALID expression like the C models, but that is too expensive. 237 I = form() 238 if np.isnan(I).any(): continue 239 240 # update value and norm 241 weight *= spherical_correction 242 total += weight * I 243 norm += weight * form_volume() 244 245 if norm > 0.0: 246 return (scale/norm)*total + background 206 247 else: 207 stride = np.array([1]) 208 vol_weight_index = slice(None, None) 209 # keep lint happy 210 fast_value = [None] 211 fast_weight = [None] 212 213 ret = np.zeros_like(args[0]) 214 norm = np.zeros_like(ret) 215 vol = np.zeros_like(ret) 216 vol_norm = np.zeros_like(ret) 217 for k in range(stride[-1]): 218 # update polydispersity parameter values 219 fast_index = k % stride[0] 220 if fast_index == 0: # bottom loop complete ... check all other loops 221 if weight.size > 0: 222 for i, index, in enumerate(k % stride): 223 args[pd_index[i]] = pd[i][0][index] 224 weight[i] = pd[i][1][index] 225 else: 226 args[pd_index[0]] = fast_value[fast_index] 227 weight[0] = fast_weight[fast_index] 228 229 # Computes the weight, and if it is not sufficient then ignore this 230 # parameter set. 231 # Note: could precompute w1*...*wn so we only need to multiply by w0 232 w = np.prod(weight) 233 if w > cutoff: 234 # Note: can precompute spherical correction if theta_index is not 235 # the fast index. Correction factor for spherical integration 236 #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0 237 spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2 238 if theta_index >= 0 else 1.0) 239 #spherical_correction = 1.0 240 241 # Call the scattering function and adds it to the total. 242 I = form(*args) 243 if np.isnan(I).any(): continue 244 ret += w * I * spherical_correction 245 norm += w 246 247 # Volume normalization. 248 # If there are "volume" polydispersity parameters, then these 249 # will be used to call the form_volume function from the user 250 # supplied kernel, and accumulate a normalized weight. 251 # Note: can precompute volume norm if fast index is not a volume 252 if form_volume: 253 vol_args = [args[index] for index in vol_index] 254 vol_weight = np.prod(weight[vol_weight_index]) 255 vol += vol_weight * form_volume(*vol_args) 256 vol_norm += vol_weight 257 258 positive = (vol * vol_norm != 0.0) 259 ret[positive] *= vol_norm[positive] / vol[positive] 260 result = scale * ret / norm + background 261 return result 248 return np.ones(nq, 'd')*background 249 250 251 def _create_default_functions(model_info): 252 """ 253 Autogenerate missing functions, such as Iqxy from Iq. 254 255 This only works for Iqxy when Iq is written in python. :func:`make_source` 256 performs a similar role for Iq written in C. This also vectorizes 257 any functions that are not already marked as vectorized. 258 """ 259 _create_vector_Iq(model_info) 260 _create_vector_Iqxy(model_info) # call create_vector_Iq() first 261 262 263 def _create_vector_Iq(model_info): 264 """ 265 Define Iq as a vector function if it exists. 266 """ 267 Iq = model_info.Iq 268 if callable(Iq) and not getattr(Iq, 'vectorized', False): 269 #print("vectorizing Iq") 270 def vector_Iq(q, *args): 271 """ 272 Vectorized 1D kernel. 273 """ 274 return np.array([Iq(qi, *args) for qi in q]) 275 vector_Iq.vectorized = True 276 model_info.Iq = vector_Iq 277 278 def _create_vector_Iqxy(model_info): 279 """ 280 Define Iqxy as a vector function if it exists, or default it from Iq(). 281 """ 282 Iq, Iqxy = model_info.Iq, model_info.Iqxy 283 if callable(Iqxy): 284 if not getattr(Iqxy, 'vectorized', False): 285 #print("vectorizing Iqxy") 286 def vector_Iqxy(qx, qy, *args): 287 """ 288 Vectorized 2D kernel. 289 """ 290 return np.array([Iqxy(qxi, qyi, *args) for qxi, qyi in zip(qx, qy)]) 291 vector_Iqxy.vectorized = True 292 model_info.Iqxy = vector_Iqxy 293 elif callable(Iq): 294 #print("defaulting Iqxy") 295 # Iq is vectorized because create_vector_Iq was already called. 296 def default_Iqxy(qx, qy, *args): 297 """ 298 Default 2D kernel. 299 """ 300 return Iq(np.sqrt(qx**2 + qy**2), *args) 301 default_Iqxy.vectorized = True 302 model_info.Iqxy = default_Iqxy 303 -
sasmodels/list_pars.py
ra84a0ca r6d6508e 25 25 for name in sorted(MODELS): 26 26 model_info = load_model_info(name) 27 for p in model_info ['parameters']:27 for p in model_info.parameters.kernel_parameters: 28 28 partable.setdefault(p.name, []) 29 29 partable[p.name].append(name) -
sasmodels/mixture.py
r72a081d r7ae2b7f 12 12 """ 13 13 from copy import copy 14 import numpy as np 14 import numpy as np # type: ignore 15 15 16 from .generate import process_parameters, COMMON_PARAMETERS, Parameter 16 from .modelinfo import Parameter, ParameterTable, ModelInfo 17 from .kernel import KernelModel, Kernel 17 18 18 SCALE=0 19 BACKGROUND=1 20 EFFECT_RADIUS=2 21 VOLFRACTION=3 19 try: 20 from typing import List 21 from .details import CallDetails 22 except ImportError: 23 pass 22 24 23 25 def make_mixture_info(parts): 26 # type: (List[ModelInfo]) -> ModelInfo 24 27 """ 25 28 Create info block for product model. … … 27 30 flatten = [] 28 31 for part in parts: 29 if part ['composition'] and part['composition'][0] == 'mixture':30 flatten.extend(part ['compostion'][1])32 if part.composition and part.composition[0] == 'mixture': 33 flatten.extend(part.composition[1]) 31 34 else: 32 35 flatten.append(part) … … 34 37 35 38 # Build new parameter list 36 pars = COMMON_PARAMETERS +[]39 combined_pars = [] 37 40 for k, part in enumerate(parts): 38 41 # Parameter prefix per model, A_, B_, ... 42 # Note that prefix must also be applied to id and length_control 43 # to support vector parameters 39 44 prefix = chr(ord('A')+k) + '_' 40 for p in part['parameters']: 41 # No background on the individual mixture elements 42 if p.name == 'background': 43 continue 44 # TODO: promote Parameter to a full class 45 # this code knows too much about the implementation! 46 p = list(p) 47 if p[0] == 'scale': # allow model subtraction 48 p[3] = [-np.inf, np.inf] # limits 49 p[0] = prefix+p[0] # relabel parameters with A_, ... 50 pars.append(p) 45 combined_pars.append(Parameter(prefix+'scale')) 46 for p in part.parameters.kernel_parameters: 47 p = copy(p) 48 p.name = prefix + p.name 49 p.id = prefix + p.id 50 if p.length_control is not None: 51 p.length_control = prefix + p.length_control 52 combined_pars.append(p) 53 parameters = ParameterTable(combined_pars) 51 54 52 model_info = {}53 model_info ['id'] = '+'.join(part['id'])54 model_info ['name'] = ' + '.join(part['name'])55 model_info ['filename']= None56 model_info ['title'] = 'Mixture model with ' + model_info['name']57 model_info ['description'] = model_info['title']58 model_info ['docs'] = model_info['title']59 model_info ['category']= "custom"60 model_info ['parameters'] = pars61 #model_info ['single']= any(part['single'] for part in parts)62 model_info ['structure_factor']= False63 model_info ['variant_info']= None64 #model_info ['tests']= []65 #model_info ['source']= []55 model_info = ModelInfo() 56 model_info.id = '+'.join(part.id for part in parts) 57 model_info.name = ' + '.join(part.name for part in parts) 58 model_info.filename = None 59 model_info.title = 'Mixture model with ' + model_info.name 60 model_info.description = model_info.title 61 model_info.docs = model_info.title 62 model_info.category = "custom" 63 model_info.parameters = parameters 64 #model_info.single = any(part['single'] for part in parts) 65 model_info.structure_factor = False 66 model_info.variant_info = None 67 #model_info.tests = [] 68 #model_info.source = [] 66 69 # Iq, Iqxy, form_volume, ER, VR and sesans 67 70 # Remember the component info blocks so we can build the model 68 model_info['composition'] = ('mixture', parts) 69 process_parameters(model_info) 70 return model_info 71 model_info.composition = ('mixture', parts) 71 72 72 73 73 class MixtureModel( object):74 class MixtureModel(KernelModel): 74 75 def __init__(self, model_info, parts): 76 # type: (ModelInfo, List[KernelModel]) -> None 75 77 self.info = model_info 76 78 self.parts = parts 77 79 78 80 def __call__(self, q_vectors): 81 # type: (List[np.ndarray]) -> MixtureKernel 79 82 # Note: may be sending the q_vectors to the n times even though they 80 83 # are only needed once. It would mess up modularity quite a bit to … … 83 86 # in opencl; or both in opencl, but one in single precision and the 84 87 # other in double precision). 85 kernels = [part (q_vectors) for part in self.parts]88 kernels = [part.make_kernel(q_vectors) for part in self.parts] 86 89 return MixtureKernel(self.info, kernels) 87 90 88 91 def release(self): 92 # type: () -> None 89 93 """ 90 94 Free resources associated with the model. … … 94 98 95 99 96 class MixtureKernel( object):100 class MixtureKernel(Kernel): 97 101 def __init__(self, model_info, kernels): 98 dim = '2d' if kernels[0].q_input.is_2d else '1d' 102 # type: (ModelInfo, List[Kernel]) -> None 103 self.dim = kernels[0].dim 104 self.info = model_info 105 self.kernels = kernels 99 106 100 # fixed offsets starts at 2 for scale and background 101 fixed_pars, pd_pars = [], [] 102 offsets = [[2, 0]] 103 #vol_index = [] 104 def accumulate(fixed, pd, volume): 105 # subtract 1 from fixed since we are removing background 106 fixed_offset, pd_offset = offsets[-1] 107 #vol_index.extend(k+pd_offset for k,v in pd if v in volume) 108 offsets.append([fixed_offset + len(fixed) - 1, pd_offset + len(pd)]) 109 pd_pars.append(pd) 110 if dim == '2d': 111 for p in kernels: 112 partype = p.info['partype'] 113 accumulate(partype['fixed-2d'], partype['pd-2d'], partype['volume']) 114 else: 115 for p in kernels: 116 partype = p.info['partype'] 117 accumulate(partype['fixed-1d'], partype['pd-1d'], partype['volume']) 118 119 #self.vol_index = vol_index 120 self.offsets = offsets 121 self.fixed_pars = fixed_pars 122 self.pd_pars = pd_pars 123 self.info = model_info 124 self.kernels = kernels 125 self.results = None 126 127 def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): 128 scale, background = fixed_pars[0:2] 107 def __call__(self, call_details, value, weight, cutoff): 108 # type: (CallDetails, np.ndarray, np.ndarry, float) -> np.ndarray 109 scale, background = value[0:2] 129 110 total = 0.0 130 self.results = [] # remember the parts for plotting later 131 for k in range(len(self.offsets)-1): 132 start_fixed, start_pd = self.offsets[k] 133 end_fixed, end_pd = self.offsets[k+1] 134 part_fixed = [fixed_pars[start_fixed], 0.0] + fixed_pars[start_fixed+1:end_fixed] 135 part_pd = [pd_pars[start_pd], 0.0] + pd_pars[start_pd+1:end_pd] 136 part_result = self.kernels[k](part_fixed, part_pd) 111 # remember the parts for plotting later 112 self.results = [] 113 for kernel, kernel_details in zip(self.kernels, call_details.parts): 114 part_result = kernel(kernel_details, value, weight, cutoff) 137 115 total += part_result 138 self.results.append( scale*sum+background)116 self.results.append(part_result) 139 117 140 118 return scale*total + background 141 119 142 120 def release(self): 143 self.p_kernel.release() 144 self.q_kernel.release() 121 # type: () -> None 122 for k in self.kernels: 123 k.release() 145 124 -
sasmodels/model_test.py
r4d76711 r7ae2b7f 43 43 Precision defaults to 5 digits (relative). 44 44 """ 45 #TODO: rename to tests so that tab completion works better for models directory 46 45 47 from __future__ import print_function 46 48 … … 48 50 import unittest 49 51 50 import numpy as np 52 import numpy as np # type: ignore 51 53 52 54 from .core import list_models, load_model_info, build_model, HAVE_OPENCL 53 from .core import call_kernel, call_ER, call_VR 55 from .details import dispersion_mesh 56 from .direct_model import call_kernel, get_weights 54 57 from .exception import annotate_exception 55 56 #TODO: rename to tests so that tab completion works better for models directory 58 from .modelinfo import expand_pars 59 60 try: 61 from typing import List, Iterator, Callable 62 except ImportError: 63 pass 64 else: 65 from .modelinfo import ParameterTable, ParameterSet, TestCondition, ModelInfo 66 from .kernelpy import PyModel, PyInput, PyKernel, DType 67 68 def call_ER(model_info, pars): 69 # type: (ModelInfo, ParameterSet) -> float 70 """ 71 Call the model ER function using *values*. 72 73 *model_info* is either *model.info* if you have a loaded model, 74 or *kernel.info* if you have a model kernel prepared for evaluation. 75 """ 76 if model_info.ER is None: 77 return 1.0 78 else: 79 value, weight = _vol_pars(model_info, pars) 80 individual_radii = model_info.ER(*value) 81 return np.sum(weight*individual_radii) / np.sum(weight) 82 83 def call_VR(model_info, pars): 84 # type: (ModelInfo, ParameterSet) -> float 85 """ 86 Call the model VR function using *pars*. 87 88 *model_info* is either *model.info* if you have a loaded model, 89 or *kernel.info* if you have a model kernel prepared for evaluation. 90 """ 91 if model_info.VR is None: 92 return 1.0 93 else: 94 value, weight = _vol_pars(model_info, pars) 95 whole, part = model_info.VR(*value) 96 return np.sum(weight*part)/np.sum(weight*whole) 97 98 def _vol_pars(model_info, pars): 99 # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray] 100 vol_pars = [get_weights(p, pars) 101 for p in model_info.parameters.call_parameters 102 if p.type == 'volume'] 103 value, weight = dispersion_mesh(model_info, vol_pars) 104 return value, weight 105 57 106 58 107 def make_suite(loaders, models): 108 # type: (List[str], List[str]) -> unittest.TestSuite 59 109 """ 60 110 Construct the pyunit test suite. … … 66 116 *models* is the list of models to test, or *["all"]* to test all models. 67 117 """ 68 69 ModelTestCase = _hide_model_case_from_nosetests() 118 ModelTestCase = _hide_model_case_from_nose() 70 119 suite = unittest.TestSuite() 71 120 … … 86 135 # don't try to call cl kernel since it will not be 87 136 # available in some environmentes. 88 is_py = callable(model_info ['Iq'])137 is_py = callable(model_info.Iq) 89 138 90 139 if is_py: # kernel implemented in python … … 124 173 125 174 126 def _hide_model_case_from_nosetests(): 175 def _hide_model_case_from_nose(): 176 # type: () -> type 127 177 class ModelTestCase(unittest.TestCase): 128 178 """ … … 135 185 def __init__(self, test_name, model_info, test_method_name, 136 186 platform, dtype): 187 # type: (str, ModelInfo, str, str, DType) -> None 137 188 self.test_name = test_name 138 189 self.info = model_info … … 140 191 self.dtype = dtype 141 192 142 setattr(self, test_method_name, self. _runTest)193 setattr(self, test_method_name, self.run_all) 143 194 unittest.TestCase.__init__(self, test_method_name) 144 195 145 def _runTest(self): 196 def run_all(self): 197 # type: () -> None 146 198 smoke_tests = [ 147 [{}, 0.1, None],148 [{}, (0.1, 0.1), None],149 [{}, 'ER', None],150 [{}, 'VR', None],199 ({}, 0.1, None), 200 ({}, (0.1, 0.1), None), 201 ({}, 'ER', None), 202 ({}, 'VR', None), 151 203 ] 152 204 153 tests = self.info ['tests']205 tests = self.info.tests 154 206 try: 155 207 model = build_model(self.info, dtype=self.dtype, 156 208 platform=self.platform) 157 209 for test in smoke_tests + tests: 158 self. _run_one_test(model, test)210 self.run_one(model, test) 159 211 160 212 if not tests and self.platform == "dll": … … 170 222 raise 171 223 172 def _run_one_test(self, model, test): 173 pars, x, y = test 224 def run_one(self, model, test): 225 # type: (PyModel, TestCondition) -> None 226 user_pars, x, y = test 227 pars = expand_pars(self.info.parameters, user_pars) 174 228 175 229 if not isinstance(y, list): … … 185 239 actual = [call_VR(model.info, pars)] 186 240 elif isinstance(x[0], tuple): 187 Qx, Qy = zip(*x)188 q_vectors = [np.array( Qx), np.array(Qy)]241 qx, qy = zip(*x) 242 q_vectors = [np.array(qx), np.array(qy)] 189 243 kernel = model.make_kernel(q_vectors) 190 244 actual = call_kernel(kernel, pars) … … 194 248 actual = call_kernel(kernel, pars) 195 249 196 self.assert Greater(len(actual),0)250 self.assertTrue(len(actual) > 0) 197 251 self.assertEqual(len(y), len(actual)) 198 252 … … 210 264 211 265 def is_near(target, actual, digits=5): 266 # type: (float, float, int) -> bool 212 267 """ 213 268 Returns true if *actual* is within *digits* significant digits of *target*. … … 218 273 219 274 def main(): 275 # type: () -> int 220 276 """ 221 277 Run tests given is sys.argv. … … 251 307 python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ... 252 308 253 If -v is included on the command line, then use verbo e output.309 If -v is included on the command line, then use verbose output. 254 310 255 311 If neither opencl nor dll is specified, then models will be tested with 256 both opencland dll; the compute target is ignored for pure python models.312 both OpenCL and dll; the compute target is ignored for pure python models. 257 313 258 314 If model1 is 'all', then all except the remaining models will be tested. … … 269 325 270 326 def model_tests(): 327 # type: () -> Iterator[Callable[[], None]] 271 328 """ 272 329 Test runner visible to nosetests. … … 276 333 tests = make_suite(['opencl', 'dll'], ['all']) 277 334 for test_i in tests: 278 yield test_i. _runTest335 yield test_i.run_all 279 336 280 337 -
sasmodels/models/cylinder.c
r26141cb r03cac08 3 3 double Iqxy(double qx, double qy, double sld, double solvent_sld, 4 4 double radius, double length, double theta, double phi); 5 6 #define INVALID(v) (v.radius<0 || v.length<0) 5 7 6 8 double form_volume(double radius, double length) … … 15 17 double length) 16 18 { 17 // TODO: return NaN if radius<0 or length<0?18 19 // precompute qr and qh to save time in the loop 19 20 const double qr = q*radius; … … 47 48 double phi) 48 49 { 49 // TODO: return NaN if radius<0 or length<0?50 50 double sn, cn; // slots to hold sincos function output 51 51 -
sasmodels/models/cylinder.py
rf247314 r7ae2b7f 82 82 """ 83 83 84 import numpy as np 85 from numpy import pi, inf 84 import numpy as np # type: ignore 85 from numpy import pi, inf # type: ignore 86 86 87 87 name = "cylinder" -
sasmodels/models/flexible_cylinder.c
re6408d0 re6408d0 1 double form_volume(double length, double kuhn_length, double radius); 2 double Iq(double q, double length, double kuhn_length, double radius, 3 double sld, double solvent_sld); 4 double Iqxy(double qx, double qy, double length, double kuhn_length, 5 double radius, double sld, double solvent_sld); 6 double flexible_cylinder_kernel(double q, double length, double kuhn_length, 7 double radius, double sld, double solvent_sld); 8 9 10 double form_volume(double length, double kuhn_length, double radius) 1 static double 2 form_volume(length, kuhn_length, radius) 11 3 { 12 4 return 1.0; 13 5 } 14 6 15 double flexible_cylinder_kernel(double q, 16 double length, 17 double kuhn_length, 18 double radius, 19 double sld, 20 double solvent_sld) 7 static double 8 Iq(double q, 9 double length, 10 double kuhn_length, 11 double radius, 12 double sld, 13 double solvent_sld) 21 14 { 22 23 const double cont = sld-solvent_sld; 24 const double qr = q*radius; 25 //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr); 26 const double crossSect = sas_J1c(qr); 27 double flex = Sk_WR(q,length,kuhn_length); 28 flex *= crossSect*crossSect; 29 flex *= M_PI*radius*radius*length; 30 flex *= cont*cont; 31 flex *= 1.0e-4; 32 return flex; 15 const double contrast = sld - solvent_sld; 16 const double cross_section = sas_J1c(q*radius); 17 const double volume = M_PI*radius*radius*length; 18 const double flex = Sk_WR(q, length, kuhn_length); 19 return 1.0e-4 * volume * square(contrast*cross_section) * flex; 33 20 } 34 35 double Iq(double q,36 double length,37 double kuhn_length,38 double radius,39 double sld,40 double solvent_sld)41 {42 43 double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);44 return result;45 }46 47 double Iqxy(double qx, double qy,48 double length,49 double kuhn_length,50 double radius,51 double sld,52 double solvent_sld)53 {54 double q;55 q = sqrt(qx*qx+qy*qy);56 double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);57 58 return result;59 } -
sasmodels/models/gel_fit.c
r30b4ddf r03cac08 1 double form_volume(void); 2 3 double Iq(double q, 4 double guinier_scale, 5 double lorentzian_scale, 6 double gyration_radius, 7 double fractal_exp, 8 double cor_length); 9 10 double Iqxy(double qx, double qy, 11 double guinier_scale, 12 double lorentzian_scale, 13 double gyration_radius, 14 double fractal_exp, 15 double cor_length); 16 17 static double _gel_fit_kernel(double q, 1 static double Iq(double q, 18 2 double guinier_scale, 19 3 double lorentzian_scale, … … 24 8 // Lorentzian Term 25 9 ////////////////////////double a(x[i]*x[i]*zeta*zeta); 26 double lorentzian_term = q*q*cor_length*cor_length;10 double lorentzian_term = square(q*cor_length); 27 11 lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term; 28 12 lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) ); … … 30 14 // Exponential Term 31 15 ////////////////////////double d(x[i]*x[i]*rg*rg); 32 double exp_term = q*q*gyration_radius*gyration_radius;16 double exp_term = square(q*gyration_radius); 33 17 exp_term = exp(-1.0*(exp_term/3.0)); 34 18 … … 37 21 return result; 38 22 } 39 double form_volume(void){40 // Unused, so free to return garbage.41 return NAN;42 }43 44 double Iq(double q,45 double guinier_scale,46 double lorentzian_scale,47 double gyration_radius,48 double fractal_exp,49 double cor_length)50 {51 return _gel_fit_kernel(q,52 guinier_scale,53 lorentzian_scale,54 gyration_radius,55 fractal_exp,56 cor_length);57 }58 59 // Iqxy is never called since no orientation or magnetic parameters.60 double Iqxy(double qx, double qy,61 double guinier_scale,62 double lorentzian_scale,63 double gyration_radius,64 double fractal_exp,65 double cor_length)66 {67 double q = sqrt(qx*qx + qy*qy);68 return _gel_fit_kernel(q,69 guinier_scale,70 lorentzian_scale,71 gyration_radius,72 fractal_exp,73 cor_length);74 }75 -
sasmodels/models/hardsphere.py
rec45c4f rd2bb604 149 149 """ 150 150 151 Iqxy = """152 // never called since no orientation or magnetic parameters.153 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);154 """155 156 151 # ER defaults to 0.0 157 152 # VR defaults to 1.0 -
sasmodels/models/hayter_msa.py
rec45c4f rd2bb604 87 87 return 1.0; 88 88 """ 89 Iqxy = """90 // never called since no orientation or magnetic parameters.91 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);92 """93 89 # ER defaults to 0.0 94 90 # VR defaults to 1.0 -
sasmodels/models/lamellar.py
rec45c4f rd2bb604 82 82 """ 83 83 84 Iqxy = """85 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);86 """87 88 84 # ER defaults to 0.0 89 85 # VR defaults to 1.0 -
sasmodels/models/lamellar_hg.py
rec45c4f rd2bb604 101 101 """ 102 102 103 Iqxy = """104 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);105 """106 107 103 # ER defaults to 0.0 108 104 # VR defaults to 1.0 -
sasmodels/models/lamellar_hg_stack_caille.py
rec45c4f rd2bb604 120 120 """ 121 121 122 Iqxy = """123 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);124 """125 126 122 # ER defaults to 0.0 127 123 # VR defaults to 1.0 -
sasmodels/models/lamellar_stack_caille.py
rec45c4f rd2bb604 104 104 """ 105 105 106 Iqxy = """107 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);108 """109 110 106 # ER defaults to 0.0 111 107 # VR defaults to 1.0 -
sasmodels/models/lamellar_stack_paracrystal.py
rec45c4f rd2bb604 132 132 """ 133 133 134 Iqxy = """135 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);136 """137 138 134 # ER defaults to 0.0 139 135 # VR defaults to 1.0 -
sasmodels/models/lib/sas_JN.c
re6408d0 re6408d0 48 48 */ 49 49 50 static double 51 sas_JN( int n, double x ) { 50 double sas_JN( int n, double x ); 51 52 double sas_JN( int n, double x ) { 52 53 53 54 const double MACHEP = 1.11022302462515654042E-16; -
sasmodels/models/lib/sph_j1c.c
re6f1410 rba32cdd 7 7 * using double precision that are the source. 8 8 */ 9 double sph_j1c(double q); 9 10 10 11 // The choice of the number of terms in the series and the cutoff value for … … 43 44 #endif 44 45 45 double sph_j1c(double q);46 46 double sph_j1c(double q) 47 47 { -
sasmodels/models/lib/sphere_form.c
rad90df9 rba32cdd 1 inline double 2 sphere_volume(double radius) 1 double sphere_volume(double radius); 2 double sphere_form(double q, double radius, double sld, double solvent_sld); 3 4 double sphere_volume(double radius) 3 5 { 4 6 return M_4PI_3*cube(radius); 5 7 } 6 8 7 inline double 8 sphere_form(double q, double radius, double sld, double solvent_sld) 9 double sphere_form(double q, double radius, double sld, double solvent_sld) 9 10 { 10 11 const double fq = sphere_volume(radius) * sph_j1c(q*radius); -
sasmodels/models/lib/wrc_cyl.c
re7678b2 rba32cdd 2 2 Functions for WRC implementation of flexible cylinders 3 3 */ 4 double Sk_WR(double q, double L, double b); 5 6 4 7 static double 5 8 AlphaSquare(double x) … … 363 366 } 364 367 365 double Sk_WR(double q, double L, double b);366 368 double Sk_WR(double q, double L, double b) 367 369 { -
sasmodels/models/onion.c
rfdb1487 rce896fd 4 4 double thickness, double A) 5 5 { 6 const double vol = 4.0/3.0 * M_PI * r * r * r;6 const double vol = M_4PI_3 * cube(r); 7 7 const double qr = q * r; 8 8 const double alpha = A * r/thickness; … … 19 19 double sinqr, cosqr; 20 20 SINCOS(qr, sinqr, cosqr); 21 fun = -3.0 (21 fun = -3.0*( 22 22 ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq 23 23 - (alpha*sinqr/qr - cosqr) … … 32 32 f_linear(double q, double r, double sld, double slope) 33 33 { 34 const double vol = 4.0/3.0 * M_PI * r * r * r;34 const double vol = M_4PI_3 * cube(r); 35 35 const double qr = q * r; 36 36 const double bes = sph_j1c(qr); … … 52 52 { 53 53 const double bes = sph_j1c(q * r); 54 const double vol = 4.0/3.0 * M_PI * r * r * r;54 const double vol = M_4PI_3 * cube(r); 55 55 return sld * vol * bes; 56 56 } … … 64 64 r += thickness[i]; 65 65 } 66 return 4.0/3.0 * M_PI * r * r * r;66 return M_4PI_3*cube(r); 67 67 } 68 68 69 69 static double 70 Iq(double q, double core_sld, double core_radius, double solvent_sld,71 double n, double in_sld[], double out_sld[], double thickness[],70 Iq(double q, double sld_core, double core_radius, double sld_solvent, 71 double n, double sld_in[], double sld_out[], double thickness[], 72 72 double A[]) 73 73 { 74 74 int i; 75 r = core_radius;76 f = f_constant(q, r, core_sld);75 double r = core_radius; 76 double f = f_constant(q, r, sld_core); 77 77 for (i=0; i<n; i++){ 78 78 const double r0 = r; … … 92 92 } 93 93 } 94 f -= f_constant(q, r, s olvent_sld);95 f2 = f * f * 1.0e-4;94 f -= f_constant(q, r, sld_solvent); 95 const double f2 = f * f * 1.0e-4; 96 96 97 97 return f2; -
sasmodels/models/onion.py
r416609b rfa5fd8d 293 293 294 294 # ["name", "units", default, [lower, upper], "type","description"], 295 parameters = [[" core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",295 parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", 296 296 "Core scattering length density"], 297 297 ["core_radius", "Ang", 200., [0, inf], "volume", 298 298 "Radius of the core"], 299 ["s olvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "",299 ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "", 300 300 "Solvent scattering length density"], 301 ["n ", "", 1, [0, 10], "volume",301 ["n_shells", "", 1, [0, 10], "volume", 302 302 "number of shells"], 303 [" in_sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",303 ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 304 304 "scattering length density at the inner radius of shell k"], 305 [" out_sld[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "",305 ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 306 306 "scattering length density at the outer radius of shell k"], 307 ["thickness[n ]", "Ang", 40., [0, inf], "volume",307 ["thickness[n_shells]", "Ang", 40., [0, inf], "volume", 308 308 "Thickness of shell k"], 309 ["A[n ]", "", 1.0, [-inf, inf], "",309 ["A[n_shells]", "", 1.0, [-inf, inf], "", 310 310 "Decay rate of shell k"], 311 311 ] 312 312 313 #source = ["lib/sph_j1c.c", "onion.c"] 314 315 def Iq(q, *args, **kw): 316 return q 317 318 def Iqxy(qx, *args, **kw): 319 return qx 320 321 322 def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A): 313 source = ["lib/sph_j1c.c", "onion.c"] 314 315 #def Iq(q, *args, **kw): 316 # return q 317 318 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 319 def profile(core_sld, core_radius, solvent_sld, n_shells, 320 in_sld, out_sld, thickness, A): 323 321 """ 324 SLD profile322 Returns shape profile with x=radius, y=SLD. 325 323 """ 326 324 327 total_radius = 1.25*(sum(thickness[:n ]) + core_radius + 1)325 total_radius = 1.25*(sum(thickness[:n_shells]) + core_radius + 1) 328 326 dr = total_radius/400 # 400 points for a smooth plot 329 327 … … 338 336 339 337 # add in the shells 340 for k in range(n ):338 for k in range(n_shells): 341 339 # Left side of each shells 342 340 r0 = r[-1] … … 365 363 beta.append(solvent_sld) 366 364 367 return np.asarray(r), np.asarray(beta) 365 return np.asarray(r), np.asarray(beta)*1e-6 368 366 369 367 def ER(core_radius, n, thickness): … … 374 372 375 373 demo = { 376 "s olvent_sld": 2.2,377 " core_sld": 1.0,374 "sld_solvent": 2.2, 375 "sld_core": 1.0, 378 376 "core_radius": 100, 379 377 "n": 4, 380 " in_sld": [0.5, 1.5, 0.9, 2.0],381 " out_sld": [nan, 0.9, 1.2, 1.6],378 "sld_in": [0.5, 1.5, 0.9, 2.0], 379 "sld_out": [nan, 0.9, 1.2, 1.6], 382 380 "thickness": [50, 75, 150, 75], 383 381 "A": [0, -1, 1e-4, 1], -
sasmodels/models/rpa.c
r13ed84c rd2bb604 1 1 double Iq(double q, double case_num, 2 double Na, double Phia, double va, double a_sld, double ba, 3 double Nb, double Phib, double vb, double b_sld, double bb, 4 double Nc, double Phic, double vc, double c_sld, double bc, 5 double Nd, double Phid, double vd, double d_sld, double bd, 2 double N[], double Phi[], double v[], double L[], double b[], 6 3 double Kab, double Kac, double Kad, 7 4 double Kbc, double Kbd, double Kcd 8 5 ); 9 6 10 double Iqxy(double qx, double qy, double case_num,11 double Na, double Phia, double va, double a_sld, double ba,12 double Nb, double Phib, double vb, double b_sld, double bb,13 double Nc, double Phic, double vc, double c_sld, double bc,14 double Nd, double Phid, double vd, double d_sld, double bd,15 double Kab, double Kac, double Kad,16 double Kbc, double Kbd, double Kcd17 );18 19 double form_volume(void);20 21 double form_volume(void)22 {23 return 1.0;24 }25 26 7 double Iq(double q, double case_num, 27 double Na, double Phia, double va, double La, double ba, 28 double Nb, double Phib, double vb, double Lb, double bb, 29 double Nc, double Phic, double vc, double Lc, double bc, 30 double Nd, double Phid, double vd, double Ld, double bd, 8 double N[], double Phi[], double v[], double L[], double b[], 31 9 double Kab, double Kac, double Kad, 32 10 double Kbc, double Kbd, double Kcd … … 36 14 #if 0 // Sasview defaults 37 15 if (icase <= 1) { 38 N a=Nb=1000.0;39 Phi a=Phib=0.0000001;16 N[0]=N[1]=1000.0; 17 Phi[0]=Phi[1]=0.0000001; 40 18 Kab=Kac=Kad=Kbc=Kbd=-0.0004; 41 L a=Lb=1e-12;42 v a=vb=100.0;43 b a=bb=5.0;19 L[0]=L[1]=1e-12; 20 v[0]=v[1]=100.0; 21 b[0]=b[1]=5.0; 44 22 } else if (icase <= 4) { 45 Phi a=0.0000001;23 Phi[0]=0.0000001; 46 24 Kab=Kac=Kad=-0.0004; 47 L a=1e-12;48 v a=100.0;49 b a=5.0;25 L[0]=1e-12; 26 v[0]=100.0; 27 b[0]=5.0; 50 28 } 51 29 #else 52 30 if (icase <= 1) { 53 N a=Nb=0.0;54 Phi a=Phib=0.0;31 N[0]=N[1]=0.0; 32 Phi[0]=Phi[1]=0.0; 55 33 Kab=Kac=Kad=Kbc=Kbd=0.0; 56 L a=Lb=Ld;57 v a=vb=vd;58 b a=bb=0.0;34 L[0]=L[1]=L[3]; 35 v[0]=v[1]=v[3]; 36 b[0]=b[1]=0.0; 59 37 } else if (icase <= 4) { 60 N a= 0.0;61 Phi a=0.0;38 N[0] = 0.0; 39 Phi[0]=0.0; 62 40 Kab=Kac=Kad=0.0; 63 L a=Ld;64 v a=vd;65 b a=0.0;41 L[0]=L[3]; 42 v[0]=v[3]; 43 b[0]=0.0; 66 44 } 67 45 #endif 68 46 69 const double Xa = q*q*b a*ba*Na/6.0;70 const double Xb = q*q*b b*bb*Nb/6.0;71 const double Xc = q*q*b c*bc*Nc/6.0;72 const double Xd = q*q*b d*bd*Nd/6.0;47 const double Xa = q*q*b[0]*b[0]*N[0]/6.0; 48 const double Xb = q*q*b[1]*b[1]*N[1]/6.0; 49 const double Xc = q*q*b[2]*b[2]*N[2]/6.0; 50 const double Xd = q*q*b[3]*b[3]*N[3]/6.0; 73 51 74 52 // limit as Xa goes to 0 is 1 … … 98 76 #if 0 99 77 const double S0aa = icase<5 100 ? 1.0 : N a*Phia*va*Paa;78 ? 1.0 : N[0]*Phi[0]*v[0]*Paa; 101 79 const double S0bb = icase<2 102 ? 1.0 : N b*Phib*vb*Pbb;103 const double S0cc = N c*Phic*vc*Pcc;104 const double S0dd = N d*Phid*vd*Pdd;80 ? 1.0 : N[1]*Phi[1]*v[1]*Pbb; 81 const double S0cc = N[2]*Phi[2]*v[2]*Pcc; 82 const double S0dd = N[3]*Phi[3]*v[3]*Pdd; 105 83 const double S0ab = icase<8 106 ? 0.0 : sqrt(N a*va*Phia*Nb*vb*Phib)*Pa*Pb;84 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 107 85 const double S0ac = icase<9 108 ? 0.0 : sqrt(N a*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);86 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 109 87 const double S0ad = icase<9 110 ? 0.0 : sqrt(N a*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);88 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 111 89 const double S0bc = (icase!=4 && icase!=7 && icase!= 9) 112 ? 0.0 : sqrt(N b*vb*Phib*Nc*vc*Phic)*Pb*Pc;90 ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 113 91 const double S0bd = (icase!=4 && icase!=7 && icase!= 9) 114 ? 0.0 : sqrt(N b*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);92 ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 115 93 const double S0cd = (icase==0 || icase==2 || icase==5) 116 ? 0.0 : sqrt(N c*vc*Phic*Nd*vd*Phid)*Pc*Pd;94 ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 117 95 #else // sasview equivalent 118 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N c,Phic,vc,Pcc);119 double S0aa = N a*Phia*va*Paa;120 double S0bb = N b*Phib*vb*Pbb;121 double S0cc = N c*Phic*vc*Pcc;122 double S0dd = N d*Phid*vd*Pdd;123 double S0ab = sqrt(N a*va*Phia*Nb*vb*Phib)*Pa*Pb;124 double S0ac = sqrt(N a*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);125 double S0ad = sqrt(N a*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);126 double S0bc = sqrt(N b*vb*Phib*Nc*vc*Phic)*Pb*Pc;127 double S0bd = sqrt(N b*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);128 double S0cd = sqrt(N c*vc*Phic*Nd*vd*Phid)*Pc*Pd;96 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc); 97 double S0aa = N[0]*Phi[0]*v[0]*Paa; 98 double S0bb = N[1]*Phi[1]*v[1]*Pbb; 99 double S0cc = N[2]*Phi[2]*v[2]*Pcc; 100 double S0dd = N[3]*Phi[3]*v[3]*Pdd; 101 double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 102 double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 103 double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 104 double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 105 double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 106 double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 129 107 switch(icase){ 130 108 case 0: … … 311 289 // Note: 1e-13 to convert from fm to cm for scattering length 312 290 const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; 313 const double Lad = icase<5 ? 0.0 : (L a/va - Ld/vd)*sqrt_Nav;314 const double Lbd = icase<2 ? 0.0 : (L b/vb - Ld/vd)*sqrt_Nav;315 const double Lcd = (L c/vc - Ld/vd)*sqrt_Nav;291 const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav; 292 const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav; 293 const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav; 316 294 317 295 const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 … … 321 299 322 300 } 323 324 double Iqxy(double qx, double qy,325 double case_num,326 double Na, double Phia, double va, double a_sld, double ba,327 double Nb, double Phib, double vb, double b_sld, double bb,328 double Nc, double Phic, double vc, double c_sld, double bc,329 double Nd, double Phid, double vd, double d_sld, double bd,330 double Kab, double Kac, double Kad,331 double Kbc, double Kbd, double Kcd332 )333 {334 double q = sqrt(qx*qx + qy*qy);335 return Iq(q,336 case_num,337 Na, Phia, va, a_sld, ba,338 Nb, Phib, vb, b_sld, bb,339 Nc, Phic, vc, c_sld, bc,340 Nd, Phid, vd, d_sld, bd,341 Kab, Kac, Kad,342 Kbc, Kbd, Kcd);343 } -
sasmodels/models/rpa.py
rec45c4f ra5b8477 86 86 # ["name", "units", default, [lower, upper], "type","description"], 87 87 parameters = [ 88 ["case_num", CASES, 0, [0, 10], "", "Component organization"],88 ["case_num", "", 1, [CASES], "", "Component organization"], 89 89 90 ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"], 91 ["Phia", "", 0.25, [0, 1], "", "volume fraction"], 92 ["va", "mL/mol", 100.0, [0, inf], "", "specific volume"], 93 ["La", "fm", 10.0, [-inf, inf], "", "scattering length"], 94 ["ba", "Ang", 5.0, [0, inf], "", "segment length"], 95 96 ["Nb", "", 1000.0, [1, inf], "", "Degree of polymerization"], 97 ["Phib", "", 0.25, [0, 1], "", "volume fraction"], 98 ["vb", "mL/mol", 100.0, [0, inf], "", "specific volume"], 99 ["Lb", "fm", 10.0, [-inf, inf], "", "scattering length"], 100 ["bb", "Ang", 5.0, [0, inf], "", "segment length"], 101 102 ["Nc", "", 1000.0, [1, inf], "", "Degree of polymerization"], 103 ["Phic", "", 0.25, [0, 1], "", "volume fraction"], 104 ["vc", "mL/mol", 100.0, [0, inf], "", "specific volume"], 105 ["Lc", "fm", 10.0, [-inf, inf], "", "scattering length"], 106 ["bc", "Ang", 5.0, [0, inf], "", "segment length"], 107 108 ["Nd", "", 1000.0, [1, inf], "", "Degree of polymerization"], 109 ["Phid", "", 0.25, [0, 1], "", "volume fraction"], 110 ["vd", "mL/mol", 100.0, [0, inf], "", "specific volume"], 111 ["Ld", "fm", 10.0, [-inf, inf], "", "scattering length"], 112 ["bd", "Ang", 5.0, [0, inf], "", "segment length"], 90 ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 91 ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 92 ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 93 ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 94 ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 113 95 114 96 ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"], -
sasmodels/models/spherical_sld.py
rec45c4f rd2bb604 170 170 # pylint: disable=bad-whitespace, line-too-long 171 171 # ["name", "units", default, [lower, upper], "type", "description"], 172 parameters = [["n _shells","", 1, [0, 9], "", "number of shells"],172 parameters = [["n", "", 1, [0, 9], "", "number of shells"], 173 173 ["radius_core", "Ang", 50.0, [0, inf], "", "intern layer thickness"], 174 174 ["sld_core", "1e-6/Ang^2", 2.07, [-inf, inf], "", "sld function flat"], … … 192 192 193 193 demo = dict( 194 n _shells=4,194 n=4, 195 195 scale=1.0, 196 196 solvent_sld=1.0, -
sasmodels/models/squarewell.py
rec45c4f rd2bb604 123 123 """ 124 124 125 Iqxy = """126 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);127 """128 129 125 # ER defaults to 0.0 130 126 # VR defaults to 1.0 -
sasmodels/models/stickyhardsphere.py
rec45c4f rd2bb604 171 171 """ 172 172 173 Iqxy = """174 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);175 """176 177 173 # ER defaults to 0.0 178 174 # VR defaults to 1.0 -
sasmodels/product.py
rf247314 r7ae2b7f 11 11 *ProductModel(P, S)*. 12 12 """ 13 import numpy as np 13 import numpy as np # type: ignore 14 14 15 from .core import call_ER_VR 16 from .generate import process_parameters 15 from .details import dispersion_mesh 16 from .modelinfo import suffix_parameter, ParameterTable, ModelInfo 17 from .kernel import KernelModel, Kernel 17 18 18 SCALE=0 19 BACKGROUND=1 20 RADIUS_EFFECTIVE=2 21 VOLFRACTION=3 19 try: 20 from typing import Tuple 21 from .modelinfo import ParameterSet 22 from .details import CallDetails 23 except ImportError: 24 pass 25 26 # TODO: make estimates available to constraints 27 #ESTIMATED_PARAMETERS = [ 28 # ["est_effect_radius", "A", 0.0, [0, np.inf], "", "Estimated effective radius"], 29 # ["est_volume_ratio", "", 1.0, [0, np.inf], "", "Estimated volume ratio"], 30 #] 22 31 23 32 # TODO: core_shell_sphere model has suppressed the volume ratio calculation 24 33 # revert it after making VR and ER available at run time as constraints. 25 34 def make_product_info(p_info, s_info): 35 # type: (ModelInfo, ModelInfo) -> ModelInfo 26 36 """ 27 37 Create info block for product model. 28 38 """ 29 p_id, p_name, p_pars = p_info ['id'], p_info['name'], p_info['parameters']30 s_id, s_name, s_pars = s_info ['id'], s_info['name'], s_info['parameters']31 # We require models to start with scale and background32 assert s_pars[SCALE].name == 'scale'33 assert s_pars[BACKGROUND].name == 'background' 34 # We require structure factors to start with effect radius and volfraction35 assert s_pars[RADIUS_EFFECTIVE].name == 'radius_effective'36 assert s_pars[VOLFRACTION].name == 'volfraction'37 # Combine the parameter sets. We are skipping the first three38 # parameters of S since scale, background are defined in P and39 # effect_radius is set from P.ER().40 pars = p_pars + s_pars[3:]41 # check for duplicates; can't use assertion since they may never be checked42 if len(set(p.name for p in pars)) != len(pars):43 raise ValueError("Duplicate parameters in %s and %s"%(p_id)) 44 model_info = {}45 model_info ['id']= '*'.join((p_id, s_id))46 model_info ['name']= ' X '.join((p_name, s_name))47 model_info ['filename']= None48 model_info ['title'] = 'Product of %s and structure factor%s'%(p_name, s_name)49 model_info ['description'] = model_info['title']50 model_info ['docs'] = model_info['title']51 model_info ['category']= "custom"52 model_info ['parameters'] = pars53 #model_info ['single'] = p_info['single'] and s_info['single']54 model_info ['structure_factor']= False55 model_info ['variant_info']= None56 #model_info ['tests']= []57 #model_info ['source']= []39 p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters 40 s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 41 p_set = set(p.id for p in p_pars.call_parameters) 42 s_set = set(p.id for p in s_pars.call_parameters) 43 44 if p_set & s_set: 45 # there is some overlap between the parameter names; tag the 46 # overlapping S parameters with name_S 47 s_list = [(suffix_parameter(par, "_S") if par.id in p_set else par) 48 for par in s_pars.kernel_parameters] 49 combined_pars = p_pars.kernel_parameters + s_list 50 else: 51 combined_pars = p_pars.kernel_parameters + s_pars.kernel_parameters 52 parameters = ParameterTable(combined_pars) 53 54 model_info = ModelInfo() 55 model_info.id = '*'.join((p_id, s_id)) 56 model_info.name = ' X '.join((p_name, s_name)) 57 model_info.filename = None 58 model_info.title = 'Product of %s and %s'%(p_name, s_name) 59 model_info.description = model_info.title 60 model_info.docs = model_info.title 61 model_info.category = "custom" 62 model_info.parameters = parameters 63 #model_info.single = p_info.single and s_info.single 64 model_info.structure_factor = False 65 model_info.variant_info = None 66 #model_info.tests = [] 67 #model_info.source = [] 58 68 # Iq, Iqxy, form_volume, ER, VR and sesans 59 model_info['composition'] = ('product', [p_info, s_info]) 60 process_parameters(model_info) 69 model_info.composition = ('product', [p_info, s_info]) 61 70 return model_info 62 71 63 class ProductModel( object):72 class ProductModel(KernelModel): 64 73 def __init__(self, model_info, P, S): 74 # type: (ModelInfo, KernelModel, KernelModel) -> None 65 75 self.info = model_info 66 76 self.P = P … … 68 78 69 79 def __call__(self, q_vectors): 80 # type: (List[np.ndarray]) -> Kernel 70 81 # Note: may be sending the q_vectors to the GPU twice even though they 71 82 # are only needed once. It would mess up modularity quite a bit to … … 74 85 # in opencl; or both in opencl, but one in single precision and the 75 86 # other in double precision). 76 p_kernel = self.P (q_vectors)77 s_kernel = self.S (q_vectors)87 p_kernel = self.P.make_kernel(q_vectors) 88 s_kernel = self.S.make_kernel(q_vectors) 78 89 return ProductKernel(self.info, p_kernel, s_kernel) 79 90 80 91 def release(self): 92 # type: (None) -> None 81 93 """ 82 94 Free resources associated with the model. … … 86 98 87 99 88 class ProductKernel( object):100 class ProductKernel(Kernel): 89 101 def __init__(self, model_info, p_kernel, s_kernel): 90 dim = '2d' if p_kernel.q_input.is_2d else '1d' 91 92 # Need to know if we want 2D and magnetic parameters when constructing 93 # a parameter map. 94 par_map = {} 95 p_info = p_kernel.info['partype'] 96 s_info = s_kernel.info['partype'] 97 vol_pars = set(p_info['volume']) 98 if dim == '2d': 99 num_p_fixed = len(p_info['fixed-2d']) 100 num_p_pd = len(p_info['pd-2d']) 101 num_s_fixed = len(s_info['fixed-2d']) 102 num_s_pd = len(s_info['pd-2d']) - 1 # exclude effect_radius 103 # volume parameters are amongst the pd pars for P, not S 104 vol_par_idx = [k for k,v in enumerate(p_info['pd-2d']) 105 if v in vol_pars] 106 else: 107 num_p_fixed = len(p_info['fixed-1d']) 108 num_p_pd = len(p_info['pd-1d']) 109 num_s_fixed = len(s_info['fixed-1d']) 110 num_s_pd = len(s_info['pd-1d']) - 1 # exclude effect_radius 111 # volume parameters are amongst the pd pars for P, not S 112 vol_par_idx = [k for k,v in enumerate(p_info['pd-1d']) 113 if v in vol_pars] 114 115 start = 0 116 par_map['p_fixed'] = np.arange(start, start+num_p_fixed) 117 # User doesn't set scale, background or effect_radius for S in P*S, 118 # so borrow values from end of p_fixed. This makes volfraction the 119 # first S parameter. 120 start += num_p_fixed 121 par_map['s_fixed'] = np.hstack(([start,start], 122 np.arange(start, start+num_s_fixed-2))) 123 par_map['volfraction'] = num_p_fixed 124 start += num_s_fixed-2 125 # vol pars offset from the start of pd pars 126 par_map['vol_pars'] = [start+k for k in vol_par_idx] 127 par_map['p_pd'] = np.arange(start, start+num_p_pd) 128 start += num_p_pd-1 129 par_map['s_pd'] = np.hstack((start, 130 np.arange(start, start+num_s_pd-1))) 131 132 self.fixed_pars = model_info['partype']['fixed-' + dim] 133 self.pd_pars = model_info['partype']['pd-' + dim] 102 # type: (ModelInfo, Kernel, Kernel) -> None 134 103 self.info = model_info 135 104 self.p_kernel = p_kernel 136 105 self.s_kernel = s_kernel 137 self.par_map = par_map138 106 139 def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): 140 pars = fixed_pars + pd_pars 141 scale = pars[SCALE] 142 background = pars[BACKGROUND] 143 s_volfraction = pars[self.par_map['volfraction']] 144 p_fixed = [pars[k] for k in self.par_map['p_fixed']] 145 s_fixed = [pars[k] for k in self.par_map['s_fixed']] 146 p_pd = [pars[k] for k in self.par_map['p_pd']] 147 s_pd = [pars[k] for k in self.par_map['s_pd']] 148 vol_pars = [pars[k] for k in self.par_map['vol_pars']] 149 107 def __call__(self, details, weights, values, cutoff): 108 # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 150 109 effect_radius, vol_ratio = call_ER_VR(self.p_kernel.info, vol_pars) 151 110 152 p_fixed[SCALE] = s_volfraction 153 p_fixed[BACKGROUND] = 0.0 154 s_fixed[SCALE] = scale 155 s_fixed[BACKGROUND] = 0.0 156 s_fixed[2] = s_volfraction/vol_ratio 157 s_pd[0] = [effect_radius], [1.0] 158 159 p_res = self.p_kernel(p_fixed, p_pd, cutoff) 160 s_res = self.s_kernel(s_fixed, s_pd, cutoff) 111 p_details, s_details = details.parts 112 p_res = self.p_kernel(p_details, weights, values, cutoff) 113 s_res = self.s_kernel(s_details, weights, values, cutoff) 161 114 #print s_fixed, s_pd, p_fixed, p_pd 162 115 163 return p_res*s_res + background116 return values[0]*(p_res*s_res) + values[1] 164 117 165 118 def release(self): 119 # type: () -> None 166 120 self.p_kernel.release() 167 self. q_kernel.release()121 self.s_kernel.release() 168 122 123 def call_ER_VR(model_info, pars): 124 """ 125 Return effect radius and volume ratio for the model. 126 """ 127 if model_info.ER is None and model_info.VR is None: 128 return 1.0, 1.0 129 130 value, weight = _vol_pars(model_info, pars) 131 132 if model_info.ER is not None: 133 individual_radii = model_info.ER(*value) 134 effect_radius = np.sum(weight*individual_radii) / np.sum(weight) 135 else: 136 effect_radius = 1.0 137 138 if model_info.VR is not None: 139 whole, part = model_info.VR(*value) 140 volume_ratio = np.sum(weight*part)/np.sum(weight*whole) 141 else: 142 volume_ratio = 1.0 143 144 return effect_radius, volume_ratio 145 146 def _vol_pars(model_info, pars): 147 # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray] 148 vol_pars = [get_weights(p, pars) 149 for p in model_info.parameters.call_parameters 150 if p.type == 'volume'] 151 value, weight = dispersion_mesh(model_info, vol_pars) 152 return value, weight 153 -
sasmodels/resolution.py
r4d76711 r7ae2b7f 6 6 from __future__ import division 7 7 8 from scipy.special import erf 9 from numpy import sqrt, log, log10 10 import numpy as np 8 from scipy.special import erf # type: ignore 9 from numpy import sqrt, log, log10, exp, pi # type: ignore 10 import numpy as np # type: ignore 11 11 12 12 __all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D", … … 35 35 smeared theory I(q). 36 36 """ 37 q = None 38 q_calc = None 37 q = None # type: np.ndarray 38 q_calc = None # type: np.ndarray 39 39 def apply(self, theory): 40 40 """ … … 476 476 *pars* are the parameter values to use when evaluating. 477 477 """ 478 from sasmodels import core478 from sasmodels import direct_model 479 479 kernel = form.make_kernel([q]) 480 theory = core.call_kernel(kernel, pars)480 theory = direct_model.call_kernel(kernel, pars) 481 481 kernel.release() 482 482 return theory … … 489 489 *q0* is the center, *dq* is the width and *q* are the points to evaluate. 490 490 """ 491 from numpy import exp, pi492 491 return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 493 492 … … 500 499 that make it slow to evaluate but give it good accuracy. 501 500 """ 502 from scipy.integrate import romberg 503 504 if any(k not in form.info['defaults'] for k in pars.keys()): 505 keys = set(form.info['defaults'].keys()) 506 extra = set(pars.keys()) - keys 507 raise ValueError("bad parameters: [%s] not in [%s]"% 508 (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 501 from scipy.integrate import romberg # type: ignore 502 503 par_set = set([p.name for p in form.info.parameters.call_parameters]) 504 if any(k not in par_set for k in pars.keys()): 505 extra = set(pars.keys()) - par_set 506 raise ValueError("bad parameters: [%s] not in [%s]" 507 % (", ".join(sorted(extra)), 508 ", ".join(sorted(pars.keys())))) 509 509 510 510 if np.isscalar(width): … … 554 554 that make it slow to evaluate but give it good accuracy. 555 555 """ 556 from scipy.integrate import romberg 557 558 if any(k not in form.info['defaults'] for k in pars.keys()): 559 keys = set(form.info['defaults'].keys()) 560 extra = set(pars.keys()) - keys 561 raise ValueError("bad parameters: [%s] not in [%s]"% 562 (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 556 from scipy.integrate import romberg # type: ignore 557 558 par_set = set([p.name for p in form.info.parameters.call_parameters]) 559 if any(k not in par_set for k in pars.keys()): 560 extra = set(pars.keys()) - par_set 561 raise ValueError("bad parameters: [%s] not in [%s]" 562 % (", ".join(sorted(extra)), 563 ", ".join(sorted(pars.keys())))) 563 564 564 565 _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) … … 692 693 693 694 def _eval_sphere(self, pars, resolution): 694 from sasmodels import core695 from sasmodels import direct_model 695 696 kernel = self.model.make_kernel([resolution.q_calc]) 696 theory = core.call_kernel(kernel, pars)697 theory = direct_model.call_kernel(kernel, pars) 697 698 result = resolution.apply(theory) 698 699 kernel.release() … … 750 751 #tol = 0.001 751 752 ## The default 3 sigma and no extra points gets 1% 752 q_calc, tol = None, 0.01 753 q_calc = None # type: np.ndarray 754 tol = 0.01 753 755 resolution = Pinhole1D(q, q_width, q_calc=q_calc) 754 756 output = self._eval_sphere(pars, resolution) 755 757 if 0: # debug plot 756 import matplotlib.pyplot as plt 758 import matplotlib.pyplot as plt # type: ignore 757 759 resolution = Perfect1D(q) 758 760 source = self._eval_sphere(pars, resolution) … … 1026 1028 """ 1027 1029 import sys 1028 import xmlrunner 1030 import xmlrunner # type: ignore 1029 1031 1030 1032 suite = unittest.TestSuite() … … 1043 1045 import sys 1044 1046 from sasmodels import core 1047 from sasmodels import direct_model 1045 1048 name = sys.argv[1] if len(sys.argv) > 1 else 'cylinder' 1046 1049 … … 1063 1066 1064 1067 kernel = model.make_kernel([resolution.q_calc]) 1065 theory = core.call_kernel(kernel, pars)1068 theory = direct_model.call_kernel(kernel, pars) 1066 1069 Iq = resolution.apply(theory) 1067 1070 … … 1073 1076 Iq_romb = romberg_pinhole_1d(resolution.q, dq, model, pars) 1074 1077 1075 import matplotlib.pyplot as plt 1078 import matplotlib.pyplot as plt # type: ignore 1076 1079 plt.loglog(resolution.q_calc, theory, label='unsmeared') 1077 1080 plt.loglog(resolution.q, Iq, label='smeared', hold=True) … … 1108 1111 Run the resolution demos. 1109 1112 """ 1110 import matplotlib.pyplot as plt 1113 import matplotlib.pyplot as plt # type: ignore 1111 1114 plt.subplot(121) 1112 1115 demo_pinhole_1d() -
sasmodels/resolution2d.py
rd6f5da6 r7ae2b7f 7 7 from __future__ import division 8 8 9 import numpy as np 10 from numpy import pi, cos, sin, sqrt 9 import numpy as np # type: ignore 10 from numpy import pi, cos, sin, sqrt # type: ignore 11 11 12 12 from . import resolution -
sasmodels/sasview_model.py
r81ec7c8 r60f03de 21 21 import logging 22 22 23 import numpy as np 23 import numpy as np # type: ignore 24 24 25 25 from . import core 26 26 from . import custom 27 27 from . import generate 28 from . import weights 29 from . import details 30 from . import modelinfo 28 31 29 32 try: 30 from typing import Dict, Mapping, Any, Sequence, Tuple, NamedTuple, List, Optional 33 from typing import Dict, Mapping, Any, Sequence, Tuple, NamedTuple, List, Optional, Union, Callable 34 from .modelinfo import ModelInfo, Parameter 31 35 from .kernel import KernelModel 32 36 MultiplicityInfoType = NamedTuple( … … 34 38 [("number", int), ("control", str), ("choices", List[str]), 35 39 ("x_axis_label", str)]) 40 SasviewModelType = Callable[[int], "SasviewModel"] 36 41 except ImportError: 37 42 pass 38 43 39 44 # TODO: separate x_axis_label from multiplicity info 40 # The x-axis label belongs with the profile generating function45 # The profile x-axis label belongs with the profile generating function 41 46 MultiplicityInfo = collections.namedtuple( 42 47 'MultiplicityInfo', … … 44 49 ) 45 50 51 # TODO: figure out how to say that the return type is a subclass 46 52 def load_standard_models(): 53 # type: () -> List[SasviewModelType] 47 54 """ 48 55 Load and return the list of predefined models. … … 55 62 try: 56 63 models.append(_make_standard_model(name)) 57 except :64 except Exception: 58 65 logging.error(traceback.format_exc()) 59 66 return models … … 61 68 62 69 def load_custom_model(path): 70 # type: (str) -> SasviewModelType 63 71 """ 64 72 Load a custom model given the model path. 65 73 """ 66 74 kernel_module = custom.load_custom_kernel_module(path) 67 model_info = generate.make_model_info(kernel_module)75 model_info = modelinfo.make_model_info(kernel_module) 68 76 return _make_model_from_info(model_info) 69 77 70 78 71 79 def _make_standard_model(name): 80 # type: (str) -> SasviewModelType 72 81 """ 73 82 Load the sasview model defined by *name*. … … 78 87 """ 79 88 kernel_module = generate.load_kernel_module(name) 80 model_info = generate.make_model_info(kernel_module)89 model_info = modelinfo.make_model_info(kernel_module) 81 90 return _make_model_from_info(model_info) 82 91 83 92 84 93 def _make_model_from_info(model_info): 94 # type: (ModelInfo) -> SasviewModelType 85 95 """ 86 96 Convert *model_info* into a SasView model wrapper. 87 97 """ 88 model_info['variant_info'] = None # temporary hack for older sasview 89 def __init__(self, multiplicity=1): 98 def __init__(self, multiplicity=None): 90 99 SasviewModel.__init__(self, multiplicity=multiplicity) 91 100 attrs = _generate_model_attributes(model_info) 92 101 attrs['__init__'] = __init__ 93 ConstructedModel = type(model_info ['id'], (SasviewModel,), attrs)102 ConstructedModel = type(model_info.name, (SasviewModel,), attrs) # type: SasviewModelType 94 103 return ConstructedModel 95 104 … … 104 113 All the attributes should be immutable to avoid accidents. 105 114 """ 115 116 # TODO: allow model to override axis labels input/output name/unit 117 118 # Process multiplicity 119 non_fittable = [] # type: List[str] 120 xlabel = model_info.profile_axes[0] if model_info.profile is not None else "" 121 variants = MultiplicityInfo(0, "", [], xlabel) 122 for p in model_info.parameters.kernel_parameters: 123 if p.name == model_info.control: 124 non_fittable.append(p.name) 125 variants = MultiplicityInfo( 126 len(p.choices), p.name, p.choices, xlabel 127 ) 128 break 129 elif p.is_control: 130 non_fittable.append(p.name) 131 variants = MultiplicityInfo( 132 int(p.limits[1]), p.name, p.choices, xlabel 133 ) 134 break 135 136 # Organize parameter sets 137 orientation_params = [] 138 magnetic_params = [] 139 fixed = [] 140 for p in model_info.parameters.user_parameters(): 141 if p.type == 'orientation': 142 orientation_params.append(p.name) 143 orientation_params.append(p.name+".width") 144 fixed.append(p.name+".width") 145 if p.type == 'magnetic': 146 orientation_params.append(p.name) 147 magnetic_params.append(p.name) 148 fixed.append(p.name+".width") 149 150 # Build class dictionary 106 151 attrs = {} # type: Dict[str, Any] 107 152 attrs['_model_info'] = model_info 108 attrs['name'] = model_info['name'] 109 attrs['id'] = model_info['id'] 110 attrs['description'] = model_info['description'] 111 attrs['category'] = model_info['category'] 112 113 # TODO: allow model to override axis labels input/output name/unit 114 115 #self.is_multifunc = False 116 non_fittable = [] # type: List[str] 117 variants = MultiplicityInfo(0, "", [], "") 118 attrs['is_structure_factor'] = model_info['structure_factor'] 119 attrs['is_form_factor'] = model_info['ER'] is not None 153 attrs['name'] = model_info.name 154 attrs['id'] = model_info.id 155 attrs['description'] = model_info.description 156 attrs['category'] = model_info.category 157 attrs['is_structure_factor'] = model_info.structure_factor 158 attrs['is_form_factor'] = model_info.ER is not None 120 159 attrs['is_multiplicity_model'] = variants[0] > 1 121 160 attrs['multiplicity_info'] = variants 122 123 partype = model_info['partype']124 orientation_params = (125 partype['orientation']126 + [n + '.width' for n in partype['orientation']]127 + partype['magnetic'])128 magnetic_params = partype['magnetic']129 fixed = [n + '.width' for n in partype['pd-2d']]130 131 161 attrs['orientation_params'] = tuple(orientation_params) 132 162 attrs['magnetic_params'] = tuple(magnetic_params) 133 163 attrs['fixed'] = tuple(fixed) 134 135 164 attrs['non_fittable'] = tuple(non_fittable) 136 165 … … 192 221 dispersion = None # type: Dict[str, Any] 193 222 #: units and limits for each parameter 194 details = None # type: Mapping[str, Tuple(str, float, float)] 195 #: multiplicity used, or None if no multiplicity controls 223 details = None # type: Dict[str, Sequence[Any]] 224 # # actual type is Dict[str, List[str, float, float]] 225 #: multiplicity value, or None if no multiplicity on the model 196 226 multiplicity = None # type: Optional[int] 197 198 def __init__(self, multiplicity): 199 # type: () -> None 200 print("initializing", self.name) 201 #raise Exception("first initialization") 202 self._model = None 203 204 ## _persistency_dict is used by sas.perspectives.fitting.basepage 205 ## to store dispersity reference. 227 #: memory for polydispersity array if using ArrayDispersion (used by sasview). 228 _persistency_dict = None # type: Dict[str, Tuple[np.ndarray, np.ndarray]] 229 230 def __init__(self, multiplicity=None): 231 # type: (Optional[int]) -> None 232 233 # TODO: _persistency_dict to persistency_dict throughout sasview 234 # TODO: refactor multiplicity to encompass variants 235 # TODO: dispersion should be a class 236 # TODO: refactor multiplicity info 237 # TODO: separate profile view from multiplicity 238 # The button label, x and y axis labels and scale need to be under 239 # the control of the model, not the fit page. Maximum flexibility, 240 # the fit page would supply the canvas and the profile could plot 241 # how it wants, but this assumes matplotlib. Next level is that 242 # we provide some sort of data description including title, labels 243 # and lines to plot. 244 245 # Get the list of hidden parameters given the mulitplicity 246 # Don't include multiplicity in the list of parameters 247 self.multiplicity = multiplicity 248 if multiplicity is not None: 249 hidden = self._model_info.get_hidden_parameters(multiplicity) 250 hidden |= set([self.multiplicity_info.control]) 251 else: 252 hidden = set() 253 206 254 self._persistency_dict = {} 207 208 self.multiplicity = multiplicity209 210 255 self.params = collections.OrderedDict() 211 256 self.dispersion = {} 212 257 self.details = {} 213 214 for p in self._model_info['parameters']: 258 for p in self._model_info.parameters.user_parameters(): 259 if p.name in hidden: 260 continue 215 261 self.params[p.name] = p.default 216 self.details[p.name] = [p.units] + p.limits 217 218 for name in self._model_info['partype']['pd-2d']: 219 self.dispersion[name] = { 220 'width': 0, 221 'npts': 35, 222 'nsigmas': 3, 223 'type': 'gaussian', 224 } 262 self.details[p.id] = [p.units, p.limits[0], p.limits[1]] 263 if p.polydisperse: 264 self.details[p.id+".width"] = [ 265 "", 0.0, 1.0 if p.relative_pd else np.inf 266 ] 267 self.dispersion[p.name] = { 268 'width': 0, 269 'npts': 35, 270 'nsigmas': 3, 271 'type': 'gaussian', 272 } 225 273 226 274 def __get_state__(self): 275 # type: () -> Dict[str, Any] 227 276 state = self.__dict__.copy() 228 277 state.pop('_model') … … 233 282 234 283 def __set_state__(self, state): 284 # type: (Dict[str, Any]) -> None 235 285 self.__dict__ = state 236 286 self._model = None 237 287 238 288 def __str__(self): 289 # type: () -> str 239 290 """ 240 291 :return: string representation … … 243 294 244 295 def is_fittable(self, par_name): 296 # type: (str) -> bool 245 297 """ 246 298 Check if a given parameter is fittable or not … … 253 305 254 306 255 # pylint: disable=no-self-use256 307 def getProfile(self): 308 # type: () -> (np.ndarray, np.ndarray) 257 309 """ 258 310 Get SLD profile … … 261 313 beta is a list of the corresponding SLD values 262 314 """ 263 return None, None 315 args = [] # type: List[Union[float, np.ndarray]] 316 for p in self._model_info.parameters.kernel_parameters: 317 if p.id == self.multiplicity_info.control: 318 args.append(float(self.multiplicity)) 319 elif p.length == 1: 320 args.append(self.params.get(p.id, np.NaN)) 321 else: 322 args.append([self.params.get(p.id+str(k), np.NaN) 323 for k in range(1,p.length+1)]) 324 return self._model_info.profile(*args) 264 325 265 326 def setParam(self, name, value): 327 # type: (str, float) -> None 266 328 """ 267 329 Set the value of a model parameter … … 290 352 291 353 def getParam(self, name): 354 # type: (str) -> float 292 355 """ 293 356 Set the value of a model parameter … … 313 376 314 377 def getParamList(self): 378 # type: () -> Sequence[str] 315 379 """ 316 380 Return a list of all available parameters for the model 317 381 """ 318 param_list = self.params.keys()382 param_list = list(self.params.keys()) 319 383 # WARNING: Extending the list with the dispersion parameters 320 384 param_list.extend(self.getDispParamList()) … … 322 386 323 387 def getDispParamList(self): 388 # type: () -> Sequence[str] 324 389 """ 325 390 Return a list of polydispersity parameters for the model 326 391 """ 327 392 # TODO: fix test so that parameter order doesn't matter 328 ret = ['%s.%s' % (d.lower(), p) 329 for d in self._model_info['partype']['pd-2d'] 330 for p in ('npts', 'nsigmas', 'width')] 393 ret = ['%s.%s' % (p.name.lower(), ext) 394 for p in self._model_info.parameters.user_parameters() 395 for ext in ('npts', 'nsigmas', 'width') 396 if p.polydisperse] 331 397 #print(ret) 332 398 return ret 333 399 334 400 def clone(self): 401 # type: () -> "SasviewModel" 335 402 """ Return a identical copy of self """ 336 403 return deepcopy(self) 337 404 338 405 def run(self, x=0.0): 406 # type: (Union[float, (float, float), List[float]]) -> float 339 407 """ 340 408 Evaluate the model … … 349 417 # pylint: disable=unpacking-non-sequence 350 418 q, phi = x 351 return self.calculate_Iq([q * math.cos(phi)], 352 [q * math.sin(phi)])[0] 353 else: 354 return self.calculate_Iq([float(x)])[0] 419 return self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)])[0] 420 else: 421 return self.calculate_Iq([x])[0] 355 422 356 423 357 424 def runXY(self, x=0.0): 425 # type: (Union[float, (float, float), List[float]]) -> float 358 426 """ 359 427 Evaluate the model in cartesian coordinates … … 366 434 """ 367 435 if isinstance(x, (list, tuple)): 368 return self.calculate_Iq([ float(x[0])], [float(x[1])])[0]369 else: 370 return self.calculate_Iq([ float(x)])[0]436 return self.calculate_Iq([x[0]], [x[1]])[0] 437 else: 438 return self.calculate_Iq([x])[0] 371 439 372 440 def evalDistribution(self, qdist): 441 # type: (Union[np.ndarray, Tuple[np.ndarray, np.ndarray], List[np.ndarray]]) -> np.ndarray 373 442 r""" 374 443 Evaluate a distribution of q-values. … … 401 470 # Check whether we have a list of ndarrays [qx,qy] 402 471 qx, qy = qdist 403 partype = self._model_info['partype'] 404 if not partype['orientation'] and not partype['magnetic']: 472 if not self._model_info.parameters.has_2d: 405 473 return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 406 474 else: … … 415 483 % type(qdist)) 416 484 417 def calculate_Iq(self, *args): 485 def calculate_Iq(self, qx, qy=None): 486 # type: (Sequence[float], Optional[Sequence[float]]) -> np.ndarray 418 487 """ 419 488 Calculate Iq for one set of q with the current parameters. … … 426 495 if self._model is None: 427 496 self._model = core.build_model(self._model_info) 428 q_vectors = [np.asarray(q) for q in args] 429 fn = self._model.make_kernel(q_vectors) 430 pars = [self.params[v] for v in fn.fixed_pars] 431 pd_pars = [self._get_weights(p) for p in fn.pd_pars] 432 result = fn(pars, pd_pars, self.cutoff) 433 fn.q_input.release() 434 fn.release() 497 if qy is not None: 498 q_vectors = [np.asarray(qx), np.asarray(qy)] 499 else: 500 q_vectors = [np.asarray(qx)] 501 kernel = self._model.make_kernel(q_vectors) 502 pairs = [self._get_weights(p) 503 for p in self._model_info.parameters.call_parameters] 504 call_details, weight, value = details.build_details(kernel, pairs) 505 result = kernel(call_details, weight, value, cutoff=self.cutoff) 506 kernel.release() 435 507 return result 436 508 437 509 def calculate_ER(self): 510 # type: () -> float 438 511 """ 439 512 Calculate the effective radius for P(q)*S(q) … … 441 514 :return: the value of the effective radius 442 515 """ 443 ER = self._model_info.get('ER', None) 444 if ER is None: 516 if self._model_info.ER is None: 445 517 return 1.0 446 518 else: 447 value s, weights= self._dispersion_mesh()448 fv = ER(*values)519 value, weight = self._dispersion_mesh() 520 fv = self._model_info.ER(*value) 449 521 #print(values[0].shape, weights.shape, fv.shape) 450 return np.sum(weight s * fv) / np.sum(weights)522 return np.sum(weight * fv) / np.sum(weight) 451 523 452 524 def calculate_VR(self): 525 # type: () -> float 453 526 """ 454 527 Calculate the volf ratio for P(q)*S(q) … … 456 529 :return: the value of the volf ratio 457 530 """ 458 VR = self._model_info.get('VR', None) 459 if VR is None: 531 if self._model_info.VR is None: 460 532 return 1.0 461 533 else: 462 value s, weights= self._dispersion_mesh()463 whole, part = VR(*values)464 return np.sum(weight s * part) / np.sum(weights* whole)534 value, weight = self._dispersion_mesh() 535 whole, part = self._model_info.VR(*value) 536 return np.sum(weight * part) / np.sum(weight * whole) 465 537 466 538 def set_dispersion(self, parameter, dispersion): 539 # type: (str, weights.Dispersion) -> Dict[str, Any] 467 540 """ 468 541 Set the dispersion object for a model parameter … … 487 560 from . import weights 488 561 disperser = weights.dispersers[dispersion.__class__.__name__] 489 dispersion = weights. models[disperser]()562 dispersion = weights.MODELS[disperser]() 490 563 self.dispersion[parameter] = dispersion.get_pars() 491 564 else: … … 493 566 494 567 def _dispersion_mesh(self): 568 # type: () -> List[Tuple[np.ndarray, np.ndarray]] 495 569 """ 496 570 Create a mesh grid of dispersion parameters and weights. … … 500 574 parameter set in the vector. 501 575 """ 502 pars = self._model_info['partype']['volume'] 503 return core.dispersion_mesh([self._get_weights(p) for p in pars]) 576 pars = [self._get_weights(p) 577 for p in self._model_info.parameters.call_parameters 578 if p.type == 'volume'] 579 return details.dispersion_mesh(self._model_info, pars) 504 580 505 581 def _get_weights(self, par): 582 # type: (Parameter) -> Tuple[np.ndarray, np.ndarray] 506 583 """ 507 584 Return dispersion weights for parameter 508 585 """ 509 from . import weights 510 relative = self._model_info['partype']['pd-rel'] 511 limits = self._model_info['limits'] 512 dis = self.dispersion[par] 513 value, weight = weights.get_weights( 514 dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 515 self.params[par], limits[par], par in relative) 516 return value, weight / np.sum(weight) 517 586 if par.name not in self.params: 587 if par.name == self.multiplicity_info.control: 588 return [self.multiplicity], [] 589 else: 590 return [np.NaN], [] 591 elif par.polydisperse: 592 dis = self.dispersion[par.name] 593 value, weight = weights.get_weights( 594 dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 595 self.params[par.name], par.limits, par.relative_pd) 596 return value, weight / np.sum(weight) 597 else: 598 return [self.params[par.name]], [] 518 599 519 600 def test_model(): 601 # type: () -> float 520 602 """ 521 603 Test that a sasview model (cylinder) can be run. … … 525 607 return cylinder.evalDistribution([0.1,0.1]) 526 608 609 def test_rpa(): 610 # type: () -> float 611 """ 612 Test that a sasview model (cylinder) can be run. 613 """ 614 RPA = _make_standard_model('rpa') 615 rpa = RPA(3) 616 return rpa.evalDistribution([0.1,0.1]) 617 527 618 528 619 def test_model_list(): 620 # type: () -> None 529 621 """ 530 622 Make sure that all models build as sasview models. -
sasmodels/sesans.py
r2684d45 r7ae2b7f 12 12 from __future__ import division 13 13 14 import numpy as np 15 from numpy import pi, exp 14 import numpy as np # type: ignore 15 from numpy import pi, exp # type: ignore 16 16 from scipy.special import jv as besselj 17 17 #import direct_model.DataMixin as model -
sasmodels/weights.py
ra936688 ra936688 3 3 """ 4 4 # TODO: include dispersion docs with the disperser models 5 from math import sqrt 6 import numpy as np 7 from scipy.special import gammaln 5 from math import sqrt # type: ignore 6 import numpy as np # type: ignore 7 from scipy.special import gammaln # type: ignore 8 8 9 9 class Dispersion(object):
Note: See TracChangeset
for help on using the changeset viewer.