# Changeset d5b5b71 in sasmodels

Ignore:
Timestamp:
Mar 1, 2017 7:29:46 AM (3 years ago)
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
b216880
Parents:
e09d1e0 (diff), e1aa129 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-814

Files:
5 edited

Unmodified
Removed
• ## sasmodels/sesans.py

 rb397165 import numpy as np  # type: ignore from numpy import pi, exp  # type: ignore from scipy.special import jv as besselj #import direct_model.DataMixin as model def make_q(q_max, Rmax): r""" Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ to $q_max$. This is the integration range of the Hankel transform; bigger range and more points makes a better numerical integration. Smaller q_min will increase reliable spin echo length range. Rmax is the "radius" of the largest expected object and can be set elsewhere. q_max is determined by the acceptance angle of the SESANS instrument. from scipy.special import j0 class SesansTransform(object): """ from sas.sascalc.data_util.nxsunit import Converter Spin-Echo SANS transform calculator.  Similar to a resolution function, the SesansTransform object takes I(q) for the set of *q_calc* values and produces a transformed dataset q_min = dq = 0.1 * 2*pi / Rmax return np.arange(q_min, Converter(q_max)(q_max, units="1/A"), dq) def make_all_q(data): *SElength* (A) is the set of spin-echo lengths in the measured data. *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin echo encoding dimension (for ToF: Q of min(R) and max(lam)). *Rmax* (A) is the maximum size sensitivity; larger radius requires more computation time. """ Return a $q$ vector suitable for calculating the total scattering cross section for calculating the effect of finite acceptance angles on Time of Flight SESANS instruments. If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed. If the instrument has a rectangular acceptance, 2 all_q vectors are needed. If the instrument has a circular acceptance, 1 all_q vector is needed """ if not data.has_no_finite_acceptance: return [] elif data.has_yz_acceptance(data): # compute qx, qy Qx, Qy = np.meshgrid(qx, qy) return [Qx, Qy] else: # else only need q # data.has_z_acceptance return [q] #: SElength from the data in the original data units; not used by transform #: but the GUI uses it, so make sure that it is present. q = None  # type: np.ndarray def transform(data, q_calc, Iq_calc, qmono, Iq_mono): """ Decides which transform type is to be used, based on the experiment data file contents (header) (2016-03-19: currently controlled from parameters script) nqmono is the number of q vectors to be used for the detector integration """ nqmono = len(qmono) if nqmono == 0: result = call_hankel(data, q_calc, Iq_calc) elif nqmono == 1: q = qmono result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) else: Qx, Qy = [qmono, qmono] Qx = np.reshape(Qx, nqx, nqy) Qy = np.reshape(Qy, nqx, nqy) Iq_mono = np.reshape(Iq_mono, nqx, nqy) qx = Qx[0, :] qy = Qy[:, 0] result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) #: q values to calculate when computing transform q_calc = None  # type: np.ndarray return result # transform arrays _H = None  # type: np.ndarray _H0 = None # type: np.ndarray def call_hankel(data, q_calc, Iq_calc): return hankel((data.x, data.x_unit), (data.lam, data.lam_unit), (data.sample.thickness, data.sample.thickness_unit), q_calc, Iq_calc) def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): return hankel(data.x, data.lam * 1e-9, data.sample.thickness / 10, q_calc, Iq_calc) def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): return hankel(data.x, data.y, data.lam * 1e-9, data.sample.thickness / 10, q_calc, Iq_calc) def TotalScatter(model, parameters):  #Work in progress!! #    Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) allq = np.linspace(0,4*pi/wavelength,1000) allIq = 1 integral = allq*allIq def __init__(self, z, SElength, zaccept, Rmax): # type: (np.ndarray, float, float) -> None #import logging; logging.info("creating SESANS transform") self.q = z self._set_hankel(SElength, zaccept, Rmax) def apply(self, Iq): # tye: (np.ndarray) -> np.ndarray G0 = np.dot(self._H0, Iq) G = np.dot(self._H.T, Iq) P = G - G0 return P def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still #============================================================================== #     2D Cosine Transform if "wavelength" is a vector #============================================================================== #allq is the q-space needed to create the total scattering cross-section def _set_hankel(self, SElength, zaccept, Rmax): # type: (np.ndarray, float, float) -> None # Force float32 arrays, otherwise run into memory problems on some machines SElength = np.asarray(SElength, dtype='float32') Gprime = np.zeros_like(wavelength, 'd') s = np.zeros_like(wavelength, 'd') sd = np.zeros_like(wavelength, 'd') Gprime = np.zeros_like(wavelength, 'd') f = np.zeros_like(wavelength, 'd') for i, wavelength_i in enumerate(wavelength): z = magfield*wavelength_i allq=np.linspace() #for calculating the Q-range of the  scattering power integral allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow alldq = (allq-allq)*1e10 sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) s[i]=1-exp(-sigma) for j, Iqy_j, qy_j in enumerate(qy): for k, Iqz_k, qz_k in enumerate(qz): Iq = np.sqrt(Iqy_j^2+Iqz_k^2) q = np.sqrt(qy_j^2 + qz_k^2) Gintegral = Iq*cos(z*Qz_k) Gprime[i] += Gintegral #                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] #                s[i] += 1-exp(Totalscatter(modelname)*thickness) #                For now, work with standard 2-phase scatter #Rmax = #value in text box somewhere in FitPage? q_max = 2*pi / (SElength - SElength) q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1]) q = np.arange(q_min, q_max, q_min, dtype='float32') dq = q_min H0 = np.float32(dq/(2*pi)) * q sd[i] += Iq f[i] = 1-s[i]+sd[i] P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] repq = np.tile(q, (SElength.size, 1)).T repSE = np.tile(SElength, (q.size, 1)) H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): #============================================================================== #     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS #============================================================================== #acceptq is the q-space needed to create limited acceptance effect SElength= wavelength*magfield G = np.zeros_like(SElength, 'd') threshold=2*pi*theta/wavelength for i, SElength_i in enumerate(SElength): allq=np.linspace() #for calculating the Q-range of the  scattering power integral allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow alldq = (allq-allq)*1e10 sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) s[i]=1-exp(-sigma) dq = (q-q)*1e10 a = (x
• ## sasmodels/compare_many.py

 r424fe00 header = ('\n"Model","%s","Count","%d","Dimension","%s"' % (name, N, "2D" if is_2d else "1D")) if not mono: header += ',"Cutoff",%g'%(cutoff,) if not mono: header += ',"Cutoff",%g'%(cutoff,) print(header) max_diff =  for k in range(N): print("%s %d"%(name, k), file=sys.stderr) print("Model %s %d"%(name, k+1), file=sys.stderr) seed = np.random.randint(1e6) pars_i = randomize_pars(model_info, pars, seed) constrain_pars(model_info, pars_i) constrain_new_to_old(model_info, pars_i) if 'sasview' in (base, comp): constrain_new_to_old(model_info, pars_i) if mono: pars_i = suppress_pd(pars_i) print("""\ MODEL is the model name of the model or "all" for all the models in alphabetical order. MODEL is the model name of the model or one of the model types listed in sasmodels.core.list_models (all, py, c, double, single, opencl, 1d, 2d, nonmagnetic, magnetic).  Model types can be combined, such as 2d+single. COUNT is the number of randomly generated parameter sets to try. A value return model = argv if not (model in MODELS) and (model != "all"): print('Bad model %s.  Use "all" or one of:'%model) target = argv try: model_list = [target] if target in MODELS else core.list_models(target) except ValueError: print('Bad model %s.  Use model type or one of:'%model) print_models() print('model types: all, py, c, double, single, opencl, 1d, 2d, nonmagnetic, magnetic') return try: data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 'accuracy': 'Low', 'view':'log', 'zero': False}) model_list = [model] if model != "all" else MODELS for model in model_list: compare_instance(model, data, index, N=count, mono=mono,
• ## sasmodels/core.py

 r52e9a45 * magnetic: models with an sld * nommagnetic: models without an sld """ if kind and kind not in KINDS: For multiple conditions, combine with plus.  For example, *c+single+2d* would return all oriented models implemented in C which can be computed accurately with single precision arithmetic. """ if kind and any(k not in KINDS for k in kind.split('+')): raise ValueError("kind not in " + ", ".join(KINDS)) files = sorted(glob(joinpath(generate.MODEL_PATH, "[a-zA-Z]*.py"))) available_models = [basename(f)[:-3] for f in files] selected = [name for name in available_models if _matches(name, kind)] if kind and '+' in kind: all_kinds = kind.split('+') condition = lambda name: all(_matches(name, k) for k in all_kinds) else: condition = lambda name: _matches(name, kind) selected = [name for name in available_models if condition(name)] return selected
• ## sasmodels/model_test.py

 r479d0f3 is_py = callable(model_info.Iq) # Some OpenCL drivers seem to be flaky, and are not producing the # expected result.  Since we don't have known test values yet for # all of our models, we are instead going to compare the results # for the 'smoke test' (that is, evaluation at q=0.1 for the default # parameters just to see that the model runs to completion) between # the OpenCL and the DLL.  To do this, we define a 'stash' which is # shared between OpenCL and DLL tests.  This is just a list.  If the # list is empty (which it will be when DLL runs, if the DLL runs # first), then the results are appended to the list.  If the list # is not empty (which it will be when OpenCL runs second), the results # are compared to the results stored in the first element of the list. # This is a horrible stateful hack which only makes sense because the # test suite is thrown away after being run once. stash = [] if is_py:  # kernel implemented in python test_name = "Model: %s, Kernel: python"%model_name test_method_name, platform="dll",  # so that dtype="double") dtype="double", stash=stash) suite.addTest(test) else:   # kernel implemented in C # test using dll if desired if 'dll' in loaders or not core.HAVE_OPENCL: test_name = "Model: %s, Kernel: dll"%model_name test_method_name = "test_%s_dll" % model_info.id test = ModelTestCase(test_name, model_info, test_method_name, platform="dll", dtype="double", stash=stash) suite.addTest(test) # test using opencl if desired and available if 'opencl' in loaders and core.HAVE_OPENCL: test = ModelTestCase(test_name, model_info, test_method_name, platform="ocl", dtype=None) platform="ocl", dtype=None, stash=stash) #print("defining", test_name) suite.addTest(test) # test using dll if desired if 'dll' in loaders or not core.HAVE_OPENCL: test_name = "Model: %s, Kernel: dll"%model_name test_method_name = "test_%s_dll" % model_info.id test = ModelTestCase(test_name, model_info, test_method_name, platform="dll", dtype="double") suite.addTest(test) """ def __init__(self, test_name, model_info, test_method_name, platform, dtype): # type: (str, ModelInfo, str, str, DType) -> None platform, dtype, stash): # type: (str, ModelInfo, str, str, DType, List[Any]) -> None self.test_name = test_name self.info = model_info self.platform = platform self.dtype = dtype self.stash = stash  # container for the results of the first run setattr(self, test_method_name, self.run_all) ] tests = self.info.tests tests = smoke_tests + self.info.tests try: model = build_model(self.info, dtype=self.dtype, platform=self.platform) for test in smoke_tests + tests: self.run_one(model, test) if not tests and self.platform == "dll": ## Uncomment the following to make forgetting the test ## values an error.  Only do so for the "dll" tests ## to reduce noise from both opencl and dll, and because ## python kernels use platform="dll". #raise Exception("No test cases provided") pass results = [self.run_one(model, test) for test in tests] if self.stash: for test, target, actual in zip(tests, self.stash, results): assert np.all(abs(target-actual)<2e-5*abs(actual)),\ "expected %s but got %s for %s"%(target, actual, test) else: self.stash.append(results) # Check for missing tests.  Only do so for the "dll" tests # to reduce noise from both opencl and dll, and because # python kernels use platform="dll". if self.platform == "dll": missing = [] ## Uncomment the following to require test cases #missing = self._find_missing_tests() if missing: raise ValueError("Missing tests for "+", ".join(missing)) except: annotate_exception(self.test_name) raise def _find_missing_tests(self): # type: () -> None """make sure there are 1D, 2D, ER and VR tests as appropriate""" model_has_VR = callable(self.info.VR) model_has_ER = callable(self.info.ER) model_has_1D = True model_has_2D = any(p.type == 'orientation' for p in self.info.parameters.kernel_parameters) # Lists of tests that have a result that is not None single = [test for test in self.info.tests if not isinstance(test, list) and test is not None] tests_has_VR = any(test == 'VR' for test in single) tests_has_ER = any(test == 'ER' for test in single) tests_has_1D_single = any(isinstance(test, float) for test in single) tests_has_2D_single = any(isinstance(test, tuple) for test in single) multiple = [test for test in self.info.tests if isinstance(test, list) and not all(result is None for result in test)] tests_has_1D_multiple = any(isinstance(test, float) for test in multiple) tests_has_2D_multiple = any(isinstance(test, tuple) for test in multiple) missing = [] if model_has_VR and not tests_has_VR: missing.append("VR") if model_has_ER and not tests_has_ER: missing.append("ER") if model_has_1D and not (tests_has_1D_single or tests_has_1D_multiple): missing.append("1D") if model_has_2D and not (tests_has_2D_single or tests_has_2D_multiple): missing.append("2D") return missing def run_one(self, model, test): if x == 'ER': actual = [call_ER(model.info, pars)] actual = np.array([call_ER(model.info, pars)]) elif x == 'VR': actual = [call_VR(model.info, pars)] actual = np.array([call_VR(model.info, pars)]) elif isinstance(x, tuple): qx, qy = zip(*x) 'f(%s); expected:%s; actual:%s' % (xi, yi, actual_yi)) return actual return ModelTestCase
• ## sasmodels/modelinfo.py

 r85fe7f8 info.docs = kernel_module.__doc__ info.category = getattr(kernel_module, 'category', None) info.single = getattr(kernel_module, 'single', True) info.opencl = getattr(kernel_module, 'opencl', True) info.structure_factor = getattr(kernel_module, 'structure_factor', False) info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) info.profile = getattr(kernel_module, 'profile', None) # type: ignore info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore # Default single and opencl to True for C models.  Python models have callable Iq. info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq)) info.single = getattr(kernel_module, 'single', not callable(info.Iq)) # multiplicity info
Note: See TracChangeset for help on using the changeset viewer.