Changeset 803f835 in sasmodels
- Timestamp:
- Jan 29, 2016 10:43:44 AM (9 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 823e620
- Parents:
- d07c883
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
extra/pylint.rc
rd4666ca r803f835 108 108 # Good variable names which should always be accepted, separated by a comma 109 109 good-names=_, 110 input, 111 h,i,j,k,n,w,x,y,z, 112 q,qx,qy,qz, 113 dt,dx,dy,dz,id, 114 Iq,dIq,Qx,Qy,Qz,Iqxy, 115 p, 110 input, id, 111 h, i, j, k, n, p, v, w, x, y, z, 112 q, qx, qy, qz, Q, Qx, Qy, Qz, 113 dt, dx, dy, dz, dq, 114 Iq, dIq, Iqxy, Iq_calc, 116 115 ER, call_ER, VR, call_VR, 117 116 Rmax, SElength, -
sasmodels/direct_model.py
r9404dd3 r803f835 1 import warnings 1 """ 2 Class interface to the model calculator. 3 4 Calling a model is somewhat non-trivial since the functions called depend 5 on the data type. For 1D data the *Iq* kernel needs to be called, for 6 2D data the *Iqxy* kernel needs to be called, and for SESANS data the 7 *Iq* kernel needs to be called followed by a Hankel transform. Before 8 the kernel is called an appropriate *q* calculation vector needs to be 9 constructed. This is not the simple *q* vector where you have measured 10 the data since the resolution calculation will require values beyond the 11 range of the measured data. After the calculation the resolution calculator 12 must be called to return the predicted value for each measured data point. 13 14 :class:`DirectModel` is a callable object that takes *parameter=value* 15 keyword arguments and returns the appropriate theory values for the data. 16 17 :class:`DataMixin` does the real work of interpreting the data and calling 18 the model calculator. This is used by :class:`DirectModel`, which uses 19 direct parameter values and by :class:`bumps_model.Experiment` which wraps 20 the parameter values in boxes so that the user can set fitting ranges, etc. 21 on the individual parameters and send the model to the Bumps optimizers. 22 """ 23 from __future__ import print_function 2 24 3 25 import numpy as np … … 19 41 currently used by :class:`sasview_model.SasviewModel` since this will 20 42 require a number of changes to SasView before we can do it. 43 44 :meth:`_interpret_data` initializes the data structures necessary 45 to manage the calculations. This sets attributes in the child class 46 such as *data_type* and *resolution*. 47 48 :meth:`_calc_theory` evaluates the model at the given control values. 49 50 :meth:`_set_data` sets the intensity data in the data object, 51 possibly with random noise added. This is useful for simulating a 52 dataset with the results from :meth:`_calc_theory`. 21 53 """ 22 54 def _interpret_data(self, data, model): 55 # pylint: disable=attribute-defined-outside-init 56 23 57 self._data = data 24 58 self._model = model … … 36 70 if self.data_type == 'sesans': 37 71 q = sesans.make_q(data.sample.zacceptance, data.Rmax) 38 self.index = slice(None, None) 72 index = slice(None, None) 73 res = None 39 74 if data.y is not None: 40 self.Iq = data.y 41 self.dIq = data.dy 75 Iq, dIq = data.y, data.dy 76 else: 77 Iq, dIq = None, None 42 78 #self._theory = np.zeros_like(q) 43 79 q_vectors = [q] … … 49 85 qmax = getattr(data, 'qmax', np.inf) 50 86 accuracy = getattr(data, 'accuracy', 'Low') 51 self.index = ~data.mask & (q >= qmin) & (q <= qmax)87 index = ~data.mask & (q >= qmin) & (q <= qmax) 52 88 if data.data is not None: 53 self.index &= ~np.isnan(data.data) 54 self.Iq = data.data[self.index] 55 self.dIq = data.err_data[self.index] 56 self.resolution = resolution2d.Pinhole2D(data=data, index=self.index, 57 nsigma=3.0, accuracy=accuracy) 89 index &= ~np.isnan(data.data) 90 Iq = data.data[index] 91 dIq = data.err_data[index] 92 else: 93 Iq, dIq = None, None 94 res = resolution2d.Pinhole2D(data=data, index=index, 95 nsigma=3.0, accuracy=accuracy) 58 96 #self._theory = np.zeros_like(self.Iq) 59 q_vectors = self.resolution.q_calc97 q_vectors = res.q_calc 60 98 elif self.data_type == 'Iq': 61 self.index = (data.x >= data.qmin) & (data.x <= data.qmax)99 index = (data.x >= data.qmin) & (data.x <= data.qmax) 62 100 if data.y is not None: 63 self.index &= ~np.isnan(data.y) 64 self.Iq = data.y[self.index] 65 self.dIq = data.dy[self.index] 101 index &= ~np.isnan(data.y) 102 Iq = data.y[index] 103 dIq = data.dy[index] 104 else: 105 Iq, dIq = None, None 66 106 if getattr(data, 'dx', None) is not None: 67 q, dq = data.x[ self.index], data.dx[self.index]68 if (dq >0).any():69 self.resolution= resolution.Pinhole1D(q, dq)107 q, dq = data.x[index], data.dx[index] 108 if (dq > 0).any(): 109 res = resolution.Pinhole1D(q, dq) 70 110 else: 71 self.resolution= resolution.Perfect1D(q)72 elif (getattr(data, 'dxl', None) is not None and73 74 self.resolution = resolution.Slit1D(data.x[self.index],75 width=data.dxh[self.index],76 height=data.dxw[self.index])77 else: 78 self.resolution = resolution.Perfect1D(data.x[self.index])111 res = resolution.Perfect1D(q) 112 elif (getattr(data, 'dxl', None) is not None 113 and getattr(data, 'dxw', None) is not None): 114 res = resolution.Slit1D(data.x[index], 115 width=data.dxh[index], 116 height=data.dxw[index]) 117 else: 118 res = resolution.Perfect1D(data.x[index]) 79 119 80 120 #self._theory = np.zeros_like(self.Iq) 81 q_vectors = [ self.resolution.q_calc]121 q_vectors = [res.q_calc] 82 122 else: 83 123 raise ValueError("Unknown data type") # never gets here … … 87 127 self._kernel_inputs = [v for v in q_vectors] 88 128 self._kernel = None 129 self.Iq, self.dIq, self.index = Iq, dIq, index 130 self.resolution = res 89 131 90 132 def _set_data(self, Iq, noise=None): 133 # pylint: disable=attribute-defined-outside-init 91 134 if noise is not None: 92 135 self.dIq = Iq*noise*0.01 … … 107 150 if self._kernel is None: 108 151 q_input = self._model.make_input(self._kernel_inputs) 109 self._kernel = self._model(q_input) 152 self._kernel = self._model(q_input) # pylint: disable=attribute-defined-outside-init 110 153 111 154 Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) … … 120 163 121 164 class DirectModel(DataMixin): 165 """ 166 Create a calculator object for a model. 167 168 *data* is 1D SAS, 2D SAS or SESANS data 169 170 *model* is a model calculator return from :func:`generate.load_model` 171 172 *cutoff* is the polydispersity weight cutoff. 173 """ 122 174 def __init__(self, data, model, cutoff=1e-5): 123 175 self.model = model 124 176 self.cutoff = cutoff 177 # Note: _interpret_data defines the model attributes 125 178 self._interpret_data(data, model) 126 179 self.kernel = make_kernel(self.model, self._kernel_inputs) 180 127 181 def __call__(self, **pars): 128 182 return self._calc_theory(pars, cutoff=self.cutoff) 183 129 184 def ER(self, **pars): 185 """ 186 Compute the equivalent radius for the model. 187 188 Return 0. if not defined. 189 """ 130 190 return call_ER(self.model.info, pars) 191 131 192 def VR(self, **pars): 193 """ 194 Compute the equivalent volume for the model, including polydispersity 195 effects. 196 197 Return 1. if not defined. 198 """ 132 199 return call_VR(self.model.info, pars) 200 133 201 def simulate_data(self, noise=None, **pars): 202 """ 203 Generate simulated data for the model. 204 """ 134 205 Iq = self.__call__(**pars) 135 206 self._set_data(Iq, noise=noise) 136 207 137 def demo(): 208 def main(): 209 """ 210 Program to evaluate a particular model at a set of q values. 211 """ 138 212 import sys 139 213 from .data import empty_data1D, empty_data2D … … 144 218 model_name = sys.argv[1] 145 219 call = sys.argv[2].upper() 146 if call not in ("ER", "VR"):220 if call not in ("ER", "VR"): 147 221 try: 148 222 values = [float(v) for v in call.split(',')] … … 153 227 data = empty_data1D([q]) 154 228 elif len(values) == 2: 155 qx, qy = values156 data = empty_data2D([qx], [qy])229 qx, qy = values 230 data = empty_data2D([qx], [qy]) 157 231 else: 158 232 print("use q or qx,qy or ER or VR") … … 164 238 model = load_model(model_definition, dtype='single') 165 239 calculator = DirectModel(data, model) 166 pars = dict((k, float(v))240 pars = dict((k, float(v)) 167 241 for pair in sys.argv[3:] 168 for k, v in [pair.split('=')])242 for k, v in [pair.split('=')]) 169 243 if call == "ER": 170 244 print(calculator.ER(**pars)) … … 176 250 177 251 if __name__ == "__main__": 178 demo()252 main()
Note: See TracChangeset
for help on using the changeset viewer.