diff -Nru sasmodels-0.98/debian/changelog sasmodels-0.99/debian/changelog --- sasmodels-0.98/debian/changelog 2018-10-06 02:30:23.000000000 +0000 +++ sasmodels-0.99/debian/changelog 2019-02-16 03:58:54.000000000 +0000 @@ -1,3 +1,17 @@ +sasmodels (0.99-2) unstable; urgency=medium + + * mark python-sasmodels-doc as Multi-Arch: foreign + + -- Drew Parsons Sat, 16 Feb 2019 11:58:54 +0800 + +sasmodels (0.99-1) experimental; urgency=medium + + * New upstream release. + * Standards-Version: 4.3.0 + * debhelper compatibility level 12 + + -- Drew Parsons Fri, 15 Feb 2019 17:53:42 +0800 + sasmodels (0.98-1) unstable; urgency=medium * New upstream release. diff -Nru sasmodels-0.98/debian/compat sasmodels-0.99/debian/compat --- sasmodels-0.98/debian/compat 2018-10-06 02:30:23.000000000 +0000 +++ sasmodels-0.99/debian/compat 1970-01-01 00:00:00.000000000 +0000 @@ -1 +0,0 @@ -11 diff -Nru sasmodels-0.98/debian/control sasmodels-0.99/debian/control --- sasmodels-0.98/debian/control 2018-10-06 02:30:23.000000000 +0000 +++ sasmodels-0.99/debian/control 2019-02-16 03:58:54.000000000 +0000 @@ -6,7 +6,7 @@ Stuart Prescott , Drew Parsons Build-Depends: - debhelper (>= 11), + debhelper-compat (= 12), dh-python, libjs-jquery, libjs-mathjax, @@ -26,7 +26,7 @@ python3-pytest, python3-scipy, python3-sphinx (>= 1.6~), -Standards-Version: 4.2.1 +Standards-Version: 4.3.0 Homepage: http://www.sasview.org/ Vcs-Git: https://salsa.debian.org/science-team/sasmodels.git Vcs-Browser: https://salsa.debian.org/science-team/sasmodels @@ -87,6 +87,7 @@ Package: python-sasmodels-doc Architecture: all +Multi-Arch: foreign Section: doc Depends: libjs-jquery, diff -Nru sasmodels-0.98/debian/patches/setup-install-so.patch sasmodels-0.99/debian/patches/setup-install-so.patch --- sasmodels-0.98/debian/patches/setup-install-so.patch 2018-10-06 02:30:23.000000000 +0000 +++ sasmodels-0.99/debian/patches/setup-install-so.patch 2019-02-16 03:58:54.000000000 +0000 @@ -3,7 +3,7 @@ --- a/setup.py +++ b/setup.py -@@ -54,11 +54,13 @@ +@@ -59,11 +59,13 @@ packages=[ 'sasmodels', 'sasmodels.models', @@ -16,5 +16,5 @@ 'sasmodels': ['*.c', '*.cl'], + 'sasmodels.compiled_models': ['*.so'], }, - install_requires=[ - ], + install_requires=install_requires, + extras_require={ diff -Nru sasmodels-0.98/doc/guide/magnetism/magnetism.rst sasmodels-0.99/doc/guide/magnetism/magnetism.rst --- sasmodels-0.98/doc/guide/magnetism/magnetism.rst 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/doc/guide/magnetism/magnetism.rst 2019-02-10 00:46:47.000000000 +0000 @@ -88,16 +88,16 @@ The user input parameters are: =========== ================================================================ - M0:sld $D_M M_0$ - mtheta:sld $\theta_M$ - mphi:sld $\phi_M$ - up:angle $\theta_\mathrm{up}$ - up:frac_i $u_i$ = (spin up)/(spin up + spin down) *before* the sample - up:frac_f $u_f$ = (spin up)/(spin up + spin down) *after* the sample + sld_M0 $D_M M_0$ + sld_mtheta $\theta_M$ + sld_mphi $\phi_M$ + up_frac_i $u_i$ = (spin up)/(spin up + spin down) *before* the sample + up_frac_f $u_f$ = (spin up)/(spin up + spin down) *after* the sample + up_angle $\theta_\mathrm{up}$ =========== ================================================================ .. note:: - The values of the 'up:frac_i' and 'up:frac_f' must be in the range 0 to 1. + The values of the 'up_frac_i' and 'up_frac_f' must be in the range 0 to 1. *Document History* diff -Nru sasmodels-0.98/doc/guide/plugin.rst sasmodels-0.99/doc/guide/plugin.rst --- sasmodels-0.98/doc/guide/plugin.rst 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/doc/guide/plugin.rst 2019-02-10 00:46:47.000000000 +0000 @@ -272,6 +272,17 @@ structure factor to account for interactions between particles. See `Form_Factors`_ for more details. +**model_info = ...** lets you define a model directly, for example, by +loading and modifying existing models. This is done implicitly by +:func:`sasmodels.core.load_model_info`, which can create a mixture model +from a pair of existing models. For example:: + + from sasmodels.core import load_model_info + model_info = load_model_info('sphere+cylinder') + +See :class:`sasmodels.modelinfo.ModelInfo` for details about the model +attributes that are defined. + Model Parameters ................ @@ -427,16 +438,16 @@ def random(): ... - -This function provides a model-specific random parameter set which shows model -features in the USANS to SANS range. For example, core-shell sphere sets the -outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q -value for the transition from flat to falling. It then uses a beta distribution -to set the percentage of the shape which is shell, giving a preference for very -thin or very thick shells (but never 0% or 100%). Using `-sets=10` in sascomp -should show a reasonable variety of curves over the default sascomp q range. -The parameter set is returned as a dictionary of `{parameter: value, ...}`. -Any model parameters not included in the dictionary will default according to + +This function provides a model-specific random parameter set which shows model +features in the USANS to SANS range. For example, core-shell sphere sets the +outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q +value for the transition from flat to falling. It then uses a beta distribution +to set the percentage of the shape which is shell, giving a preference for very +thin or very thick shells (but never 0% or 100%). Using `-sets=10` in sascomp +should show a reasonable variety of curves over the default sascomp q range. +The parameter set is returned as a dictionary of `{parameter: value, ...}`. +Any model parameters not included in the dictionary will default according to the code in the `_randomize_one()` function from sasmodels/compare.py. Python Models @@ -700,8 +711,8 @@ to test for finite and not NaN. erf, erfc, tgamma, lgamma: **do not use** Special functions that should be part of the standard, but are missing - or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma - instead (see below). Note: lgamma(x) has not yet been tested. + or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma + and sas_lgamma instead (see below). Some non-standard constants and functions are also provided: @@ -768,12 +779,30 @@ sas_gamma(x): Gamma function sas_gamma\ $(x) = \Gamma(x)$. - The standard math function, tgamma(x) is unstable for $x < 1$ + The standard math function, tgamma(x), is unstable for $x < 1$ on some platforms. :code:`source = ["lib/sas_gamma.c", ...]` (`sas_gamma.c `_) + sas_gammaln(x): + log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. + + The standard math function, lgamma(x), is incorrect for single + precision on some platforms. + + :code:`source = ["lib/sas_gammainc.c", ...]` + (`sas_gammainc.c `_) + + sas_gammainc(a, x), sas_gammaincc(a, x): + Incomplete gamma function + sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ + and complementary incomplete gamma function + sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ + + :code:`source = ["lib/sas_gammainc.c", ...]` + (`sas_gammainc.c `_) + sas_erf(x), sas_erfc(x): Error function sas_erf\ $(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$ @@ -810,6 +839,8 @@ $J_n(x) = \frac{1}{\pi}\int_0^\pi \cos(n\tau - x\sin(\tau))\,d\tau$. If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively. + Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100]. + The standard math function jn(n, x) is not available on all platforms. :code:`source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c", ...]` @@ -818,9 +849,11 @@ sas_Si(x): Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$. + Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100]. + This function uses Taylor series for small and large arguments: - For large arguments, + For large arguments use the following Taylor series, .. math:: diff -Nru sasmodels-0.98/example/multiscatfit.py sasmodels-0.99/example/multiscatfit.py --- sasmodels-0.98/example/multiscatfit.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/example/multiscatfit.py 2019-02-10 00:46:47.000000000 +0000 @@ -14,10 +14,10 @@ On Unix/Mac running as developer I do:: # Show the model without fitting - PYTHONPATH=..:../explore:../../bumps:../../sasview/src python multiscatfit.py + PYTHONPATH=..:../../bumps:../../sasview/src python multiscatfit.py # Run the fit - PYTHONPATH=..:../explore:../../bumps:../../sasview/src ../../bumps/run.py \ + PYTHONPATH=..:../../bumps:../../sasview/src ../../bumps/run.py \ multiscatfit.py --store=/tmp/t1 You may be able to run multiscatfit.py against the distributed sasview @@ -54,6 +54,12 @@ radius_equatorial_pd=.000128, radius_equatorial_pd_n=1, radius_equatorial_pd_nsigma=0, ) +# Tie the model to the data +M = Experiment(data=data, model=model) + +# Stack mulitple scattering on top of the existing resolution function. +M.resolution = MultipleScattering(resolution=M.resolution, probability=0.) + # SET THE FITTING PARAMETERS model.radius_polar.range(15, 3000) model.radius_equatorial.range(15, 3000) @@ -64,24 +70,11 @@ model.background.range(0,1000) model.scale.range(0, 0.1) -# Mulitple scattering probability parameter -# HACK: the probability is stuffed in as an extra parameter to the experiment. -probability = Parameter(name="probability", value=0.0) -probability.range(0.0, 0.9) +# The multiple scattering probability parameter is in the resolution function +# instead of the scattering function, so access it through M.resolution +M.scattering_probability.range(0.0, 0.9) -M = Experiment(data=data, model=model, extra_pars={'probability': probability}) - -# Stack mulitple scattering on top of the existing resolution function. -# Because resolution functions in sasview don't have fitting parameters, -# we instead allow the multiple scattering calculator to take a function -# instead of a probability. This function returns the current value of -# the parameter. ** THIS IS TEMPORARY ** when multiple scattering is -# properly integrated into sasmodels and sasview, its fittable parameter -# will be treated like the model parameters. -M.resolution = MultipleScattering(resolution=M.resolution, - probability=lambda: probability.value, - ) -M._kernel_inputs = M.resolution.q_calc +# Let bumps know that we are fitting this experiment problem = FitProblem(M) if __name__ == "__main__": diff -Nru sasmodels-0.98/explore/precision.py sasmodels-0.99/explore/precision.py --- sasmodels-0.98/explore/precision.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/explore/precision.py 2019-02-10 00:46:47.000000000 +0000 @@ -94,6 +94,9 @@ zoom: [1000,1010] neg: [-100,100] + For arbitrary range use "start:stop:steps:scale" where scale is + one of log, lin, or linear. + *diff* is "relative", "absolute" or "none" *x_bits* is the precision with which the x values are specified. The @@ -101,16 +104,23 @@ """ linear = not xrange.startswith("log") if xrange == "zoom": - lin_min, lin_max, lin_steps = 1000, 1010, 2000 + start, stop, steps = 1000, 1010, 2000 elif xrange == "neg": - lin_min, lin_max, lin_steps = -100.1, 100.1, 2000 + start, stop, steps = -100.1, 100.1, 2000 elif xrange == "linear": - lin_min, lin_max, lin_steps = 1, 1000, 2000 - lin_min, lin_max, lin_steps = 0.001, 2, 2000 + start, stop, steps = 1, 1000, 2000 + start, stop, steps = 0.001, 2, 2000 elif xrange == "log": - log_min, log_max, log_steps = -3, 5, 400 + start, stop, steps = -3, 5, 400 elif xrange == "logq": - log_min, log_max, log_steps = -4, 1, 400 + start, stop, steps = -4, 1, 400 + elif ':' in xrange: + parts = xrange.split(':') + linear = parts[3] != "log" if len(parts) == 4 else True + steps = int(parts[2]) if len(parts) > 2 else 400 + start = float(parts[0]) + stop = float(parts[1]) + else: raise ValueError("unknown range "+xrange) with mp.workprec(500): @@ -120,19 +130,19 @@ # rather than x to f(nearest(x)) where nearest(x) is the nearest # value to x in the given precision. if linear: - lin_min = max(lin_min, self.limits[0]) - lin_max = min(lin_max, self.limits[1]) - qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single') - #qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='double') + start = max(start, self.limits[0]) + stop = min(stop, self.limits[1]) + qrf = np.linspace(start, stop, steps, dtype='single') + #qrf = np.linspace(start, stop, steps, dtype='double') qr = [mp.mpf(float(v)) for v in qrf] - #qr = mp.linspace(lin_min, lin_max, lin_steps) + #qr = mp.linspace(start, stop, steps) else: - log_min = np.log10(max(10**log_min, self.limits[0])) - log_max = np.log10(min(10**log_max, self.limits[1])) - qrf = np.logspace(log_min, log_max, log_steps, dtype='single') - #qrf = np.logspace(log_min, log_max, log_steps, dtype='double') + start = np.log10(max(10**start, self.limits[0])) + stop = np.log10(min(10**stop, self.limits[1])) + qrf = np.logspace(start, stop, steps, dtype='single') + #qrf = np.logspace(start, stop, steps, dtype='double') qr = [mp.mpf(float(v)) for v in qrf] - #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)] + #qr = [10**v for v in mp.linspace(start, stop, steps)] target = self.call_mpmath(qr, bits=500) pylab.subplot(121) @@ -175,7 +185,7 @@ Use relative error if SHOW_DIFF, otherwise just plot the value directly. """ if diff == "relative": - err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd') + err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd') #err = np.clip(err, 0, 1) pylab.loglog(x, err, '-', label=label) elif diff == "absolute": @@ -196,6 +206,23 @@ model_info = modelinfo.make_model_info(Kernel) return model_info +# Hack to allow second parameter A in two parameter functions +A = 1 +def parse_extra_pars(): + global A + + A_str = str(A) + pop = [] + for k, v in enumerate(sys.argv[1:]): + if v.startswith("A="): + A_str = v[2:] + pop.append(k+1) + if pop: + sys.argv = [v for k, v in enumerate(sys.argv) if k not in pop] + A = float(A_str) + +parse_extra_pars() + # =============== FUNCTION DEFINITIONS ================ @@ -298,6 +325,25 @@ limits=(-3.1, 10), ) add_function( + name="gammaln(x)", + mp_function=mp.loggamma, + np_function=scipy.special.gammaln, + ocl_function=make_ocl("return sas_gammaln(q);", "sas_gammaln", ["lib/sas_gammainc.c"]), + #ocl_function=make_ocl("return lgamma(q);", "sas_gammaln"), +) +add_function( + name="gammainc(x)", + mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a), + np_function=lambda x, a=A: scipy.special.gammainc(a, x), + ocl_function=make_ocl("return sas_gammainc(%.15g,q);"%A, "sas_gammainc", ["lib/sas_gammainc.c"]), +) +add_function( + name="gammaincc(x)", + mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a), + np_function=lambda x, a=A: scipy.special.gammaincc(a, x), + ocl_function=make_ocl("return sas_gammaincc(%.15g,q);"%A, "sas_gammaincc", ["lib/sas_gammainc.c"]), +) +add_function( name="erf(x)", mp_function=mp.erf, np_function=scipy.special.erf, @@ -462,8 +508,8 @@ lanczos_gamma = """\ const double coeff[] = { - 76.18009172947146, -86.50532032941677, - 24.01409824083091, -1.231739572450155, + 76.18009172947146, -86.50532032941677, + 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2,-0.5395239384953e-5 }; const double x = q; @@ -474,7 +520,7 @@ return -tmp + log(2.5066282746310005*ser/x); """ add_function( - name="log gamma(x)", + name="loggamma(x)", mp_function=mp.loggamma, np_function=scipy.special.gammaln, ocl_function=make_ocl(lanczos_gamma, "lgamma"), @@ -598,7 +644,7 @@ ) ALL_FUNCTIONS = set(FUNCTIONS.keys()) -ALL_FUNCTIONS.discard("loggamma") # OCL version not ready yet +ALL_FUNCTIONS.discard("loggamma") # use cephes-based gammaln instead ALL_FUNCTIONS.discard("3j1/x:taylor") ALL_FUNCTIONS.discard("3j1/x:trig") ALL_FUNCTIONS.discard("2J1/x:alt") @@ -614,12 +660,17 @@ -a indicates that the absolute error should be plotted, -r indicates that the relative error should be plotted (default), -x indicates the steps in x, where is one of the following - log indicates log stepping in [10^-3, 10^5] (default) - logq indicates log stepping in [10^-4, 10^1] - linear indicates linear stepping in [1, 1000] - zoom indicates linear stepping in [1000, 1010] - neg indicates linear stepping in [-100.1, 100.1] -and name is "all" or one of: + log indicates log stepping in [10^-3, 10^5] (default) + logq indicates log stepping in [10^-4, 10^1] + linear indicates linear stepping in [1, 1000] + zoom indicates linear stepping in [1000, 1010] + neg indicates linear stepping in [-100.1, 100.1] + start:stop:n[:stepping] indicates an n-step plot in [start, stop] + or [10^start, 10^stop] if stepping is "log" (default n=400) +Some functions (notably gammainc/gammaincc) have an additional parameter A +which can be set from the command line as A=value. Default is A=1. + +Name is one of: """+names) sys.exit(1) diff -Nru sasmodels-0.98/sasmodels/bumps_model.py sasmodels-0.99/sasmodels/bumps_model.py --- sasmodels-0.98/sasmodels/bumps_model.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/bumps_model.py 2019-02-10 00:46:47.000000000 +0000 @@ -34,6 +34,7 @@ # Optional import. This allows the doc builder and nosetests to run even # when bumps is not on the path. from bumps.names import Parameter # type: ignore + from bumps.parameter import Reference # type: ignore except ImportError: pass @@ -138,12 +139,27 @@ _cache = None # type: Dict[str, np.ndarray] def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None): # type: (Data, Model, float) -> None + # Allow resolution function to define fittable parameters. We do this + # by creating reference parameters within the resolution object rather + # than modifying the object itself to use bumps parameters. We need + # to reset the parameters each time the object has changed. These + # additional parameters need to be returned from the fitting engine. + # To make them available to the user, they are added as top-level + # attributes to the experiment object. The only change to the + # resolution function is that it needs an optional 'fittable' attribute + # which maps the internal name to the user visible name for the + # for the parameter. + self._resolution = None + self._resolution_pars = {} # remember inputs so we can inspect from outside self.name = data.filename if name is None else name self.model = model self.cutoff = cutoff self._interpret_data(data, model.sasmodel) self._cache = {} + # CRUFT: no longer need extra parameters + # Multiple scattering probability is now retrieved directly from the + # multiple scattering resolution function. self.extra_pars = extra_pars def update(self): @@ -161,14 +177,38 @@ """ return len(self.Iq) + @property + def resolution(self): + return self._resolution + + @resolution.setter + def resolution(self, value): + self._resolution = value + + # Remove old resolution fitting parameters from experiment + for name in self._resolution_pars: + delattr(self, name) + + # Create new resolution fitting parameters + res_pars = getattr(self._resolution, 'fittable', {}) + self._resolution_pars = { + name: Reference(self._resolution, refname, name=name) + for refname, name in res_pars.items() + } + + # Add new resolution fitting parameters as experiment attributes + for name, ref in self._resolution_pars.items(): + setattr(self, name, ref) + def parameters(self): # type: () -> Dict[str, Parameter] """ Return a dictionary of parameters """ pars = self.model.parameters() - if self.extra_pars: + if self.extra_pars is not None: pars.update(self.extra_pars) + pars.update(self._resolution_pars) return pars def theory(self): diff -Nru sasmodels-0.98/sasmodels/compare.py sasmodels-0.99/sasmodels/compare.py --- sasmodels-0.98/sasmodels/compare.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/compare.py 2019-02-10 00:46:47.000000000 +0000 @@ -367,7 +367,7 @@ return np.random.uniform(-0.5, 12) # Limit magnetic SLDs to a smaller range, from zero to iron=5/A^2 - if par.name.startswith('M0:'): + if par.name.endswith('_M0'): return np.random.uniform(0, 5) # Guess at the random length/radius/thickness. In practice, all models @@ -537,7 +537,7 @@ magnetic = False magnetic_pars = [] for p in parameters.user_parameters(pars, is2d): - if any(p.id.startswith(x) for x in ('M0:', 'mtheta:', 'mphi:')): + if any(p.id.endswith(x) for x in ('_M0', '_mtheta', '_mphi')): continue if p.id.startswith('up:'): magnetic_pars.append("%s=%s"%(p.id, pars.get(p.id, p.default))) @@ -549,9 +549,9 @@ nsigma=pars.get(p.id+"_pd_nsgima", 3.), pdtype=pars.get(p.id+"_pd_type", 'gaussian'), relative_pd=p.relative_pd, - M0=pars.get('M0:'+p.id, 0.), - mphi=pars.get('mphi:'+p.id, 0.), - mtheta=pars.get('mtheta:'+p.id, 0.), + M0=pars.get(p.id+'_M0', 0.), + mphi=pars.get(p.id+'_mphi', 0.), + mtheta=pars.get(p.id+'_mtheta', 0.), ) lines.append(_format_par(p.name, **fields)) magnetic = magnetic or fields['M0'] != 0. @@ -617,13 +617,13 @@ pars = pars.copy() if suppress: for p in pars: - if p.startswith("M0:"): + if p.endswith("_M0"): pars[p] = 0 else: any_mag = False first_mag = None for p in pars: - if p.startswith("M0:"): + if p.endswith("_M0"): any_mag |= (pars[p] != 0) if first_mag is None: first_mag = p diff -Nru sasmodels-0.98/sasmodels/convert.py sasmodels-0.99/sasmodels/convert.py --- sasmodels-0.98/sasmodels/convert.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/convert.py 2019-02-10 00:46:47.000000000 +0000 @@ -164,8 +164,26 @@ def _hand_convert(name, oldpars, version=(3, 1, 2)): if version == (3, 1, 2): oldpars = _hand_convert_3_1_2_to_4_1(name, oldpars) + if version < (4, 2, 0): + oldpars = _rename_magnetic_pars(oldpars) return oldpars +def _rename_magnetic_pars(pars): + """ + Change from M0:par to par_M0, etc. + """ + keys = list(pars.items()) + for k in keys: + if k.startswith('M0:'): + pars[k[3:]+'_M0'] = pars.pop(k) + elif k.startswith('mtheta:'): + pars[k[7:]+'_mtheta'] = pars.pop(k) + elif k.startswith('mphi:'): + pars[k[5:]+'_mphi'] = pars.pop(k) + elif k.startswith('up:'): + pars['up_'+k[3:]] = pars.pop(k) + return pars + def _hand_convert_3_1_2_to_4_1(name, oldpars): if name == 'core_shell_parallelepiped': # Make sure pd on rim parameters defaults to zero diff -Nru sasmodels-0.98/sasmodels/direct_model.py sasmodels-0.99/sasmodels/direct_model.py --- sasmodels-0.98/sasmodels/direct_model.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/direct_model.py 2019-02-10 00:46:47.000000000 +0000 @@ -241,8 +241,6 @@ Iq, dIq = data.y, data.dy else: Iq, dIq = None, None - #self._theory = np.zeros_like(q) - q_vectors = [res.q_calc] elif self.data_type == 'Iqxy': #if not model.info.parameters.has_2d: # raise ValueError("not 2D without orientation or magnetic parameters") @@ -259,8 +257,6 @@ Iq, dIq = None, None res = resolution2d.Pinhole2D(data=data, index=index, nsigma=3.0, accuracy=accuracy) - #self._theory = np.zeros_like(self.Iq) - q_vectors = res.q_calc elif self.data_type == 'Iq': index = (data.x >= data.qmin) & (data.x <= data.qmax) mask = getattr(data, 'mask', None) @@ -285,9 +281,6 @@ qy_width=data.dxw[index]) else: res = resolution.Perfect1D(data.x[index]) - - #self._theory = np.zeros_like(self.Iq) - q_vectors = [res.q_calc] elif self.data_type == 'Iq-oriented': index = (data.x >= data.qmin) & (data.x <= data.qmax) if data.y is not None: @@ -303,13 +296,11 @@ res = resolution2d.Slit2D(data.x[index], qx_width=data.dxw[index], qy_width=data.dxl[index]) - q_vectors = res.q_calc else: raise ValueError("Unknown data type") # never gets here # Remember function inputs so we can delay loading the function and # so we can save/restore state - self._kernel_inputs = q_vectors self._kernel = None self.Iq, self.dIq, self.index = Iq, dIq, index self.resolution = res @@ -346,7 +337,13 @@ def _calc_theory(self, pars, cutoff=0.0): # type: (ParameterSet, float) -> np.ndarray if self._kernel is None: - self._kernel = self._model.make_kernel(self._kernel_inputs) + # TODO: change interfaces so that resolution returns kernel inputs + # Maybe have resolution always return a tuple, or maybe have + # make_kernel accept either an ndarray or a pair of ndarrays. + kernel_inputs = self.resolution.q_calc + if isinstance(kernel_inputs, np.ndarray): + kernel_inputs = (kernel_inputs,) + self._kernel = self._model.make_kernel(kernel_inputs) # Need to pull background out of resolution for multiple scattering background = pars.get('background', DEFAULT_BACKGROUND) diff -Nru sasmodels-0.98/sasmodels/__init__.py sasmodels-0.99/sasmodels/__init__.py --- sasmodels-0.98/sasmodels/__init__.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/__init__.py 2019-02-10 00:46:47.000000000 +0000 @@ -13,7 +13,7 @@ OpenCL drivers are available. See :mod:`generate` for details on defining new models. """ -__version__ = "0.98" +__version__ = "0.99" def data_files(): """ diff -Nru sasmodels-0.98/sasmodels/kernelpy.py sasmodels-0.99/sasmodels/kernelpy.py --- sasmodels-0.98/sasmodels/kernelpy.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/kernelpy.py 2019-02-10 00:46:47.000000000 +0000 @@ -36,6 +36,7 @@ _create_default_functions(model_info) self.info = model_info self.dtype = np.dtype('d') + logger.info("make python model " + self.info.name) def make_kernel(self, q_vectors): q_input = PyInput(q_vectors, dtype=F64) diff -Nru sasmodels-0.98/sasmodels/modelinfo.py sasmodels-0.99/sasmodels/modelinfo.py --- sasmodels-0.98/sasmodels/modelinfo.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/modelinfo.py 2019-02-10 00:46:47.000000000 +0000 @@ -465,7 +465,7 @@ for p in self.kernel_parameters) self.is_asymmetric = any(p.name == 'psi' for p in self.kernel_parameters) self.magnetism_index = [k for k, p in enumerate(self.call_parameters) - if p.id.startswith('M0:')] + if p.id.endswith('_M0')] self.pd_1d = set(p.name for p in self.call_parameters if p.polydisperse and p.type not in ('orientation', 'magnetic')) @@ -585,21 +585,21 @@ # Add the magnetic parameters to the end of the call parameter list. if self.nmagnetic > 0: full_list.extend([ - Parameter('up:frac_i', '', 0., [0., 1.], + Parameter('up_frac_i', '', 0., [0., 1.], 'magnetic', 'fraction of spin up incident'), - Parameter('up:frac_f', '', 0., [0., 1.], + Parameter('up_frac_f', '', 0., [0., 1.], 'magnetic', 'fraction of spin up final'), - Parameter('up:angle', 'degrees', 0., [0., 360.], + Parameter('up_angle', 'degrees', 0., [0., 360.], 'magnetic', 'spin up angle'), ]) slds = [p for p in full_list if p.type == 'sld'] for p in slds: full_list.extend([ - Parameter('M0:'+p.id, '1e-6/Ang^2', 0., [-np.inf, np.inf], + Parameter(p.id+'_M0', '1e-6/Ang^2', 0., [-np.inf, np.inf], 'magnetic', 'magnetic amplitude for '+p.description), - Parameter('mtheta:'+p.id, 'degrees', 0., [-90., 90.], + Parameter(p.id+'_mtheta', 'degrees', 0., [-90., 90.], 'magnetic', 'magnetic latitude for '+p.description), - Parameter('mphi:'+p.id, 'degrees', 0., [-180., 180.], + Parameter(p.id+'_mphi', 'degrees', 0., [-180., 180.], 'magnetic', 'magnetic longitude for '+p.description), ]) #print("call parameters", full_list) @@ -639,8 +639,8 @@ early, and rerender the table when it is changed. Parameters marked as sld will automatically have a set of associated - magnetic parameters (m0:p, mtheta:p, mphi:p), as well as polarization - information (up:theta, up:frac_i, up:frac_f). + magnetic parameters (p_M0, p_mtheta, p_mphi), as well as polarization + information (up_theta, up_frac_i, up_frac_f). """ # control parameters go first control = [p for p in self.kernel_parameters if p.is_control] @@ -667,9 +667,9 @@ """add the named parameter, and related magnetic parameters if any""" result.append(expanded_pars[name]) if is2d: - for tag in 'M0:', 'mtheta:', 'mphi:': - if tag+name in expanded_pars: - result.append(expanded_pars[tag+name]) + for tag in '_M0', '_mtheta', '_mphi': + if name+tag in expanded_pars: + result.append(expanded_pars[name+tag]) # Gather the user parameters in order result = control + self.COMMON @@ -702,11 +702,11 @@ else: append_group(p.id) - if is2d and 'up:angle' in expanded_pars: + if is2d and 'up_angle' in expanded_pars: result.extend([ - expanded_pars['up:frac_i'], - expanded_pars['up:frac_f'], - expanded_pars['up:angle'], + expanded_pars['up_frac_i'], + expanded_pars['up_frac_f'], + expanded_pars['up_angle'], ]) return result @@ -1015,4 +1015,9 @@ for p in self.parameters.kernel_parameters for k in range(control+1, p.length+1) if p.length > 1) + for p in self.parameters.kernel_parameters: + if p.length > 1 and p.type == "sld": + for k in range(control+1, p.length+1): + base = p.id+str(k) + hidden.update((base+"_M0", base+"_mtheta", base+"_mphi")) return hidden diff -Nru sasmodels-0.98/sasmodels/models/lib/sas_gammainc.c sasmodels-0.99/sasmodels/models/lib/sas_gammainc.c --- sasmodels-0.98/sasmodels/models/lib/sas_gammainc.c 1970-01-01 00:00:00.000000000 +0000 +++ sasmodels-0.99/sasmodels/models/lib/sas_gammainc.c 2019-02-10 00:46:47.000000000 +0000 @@ -0,0 +1,458 @@ +#if FLOAT_SIZE>4 // double precision +// based on cephes/double/igam.c +double gammaln(double x); +double gammainc(double a, double x); +double gammaincc(double a, double x); + +double cephes_igamc(double a, double x); +double cephes_igam(double a, double x); +double cephes_lgam(double x); +double cephes_lgam2(double x); + +double gammainc(double a, double x) +{ + if ((x <= 0) || ( a <= 0)) return 0.0; + if ((x > 1.0) && (x > a)) return 1.0 - cephes_igamc(a, x); + return cephes_igam(a, x); +} + +double gammaincc(double a, double x) +{ + if ((x <= 0) || (a <= 0)) return 1.0; + if ((x < 1.0) || (x < a)) return 1.0 - cephes_igam(a, x); + return cephes_igamc(a, x); +} + +double gammaln(double x) +{ + if (isnan(x)) return(x); + if (!isfinite(x)) return(INFINITY); + return cephes_lgam(x); +} + +double cephes_igamc(double a, double x) +{ + const double MACHEP = 1.11022302462515654042E-16; // IEEE 2**-53 + const double MAXLOG = 7.09782712893383996843E2; // IEEE log(2**1024) denormalized + const double BIG = 4.503599627370496e15; + const double BIGINV = 2.22044604925031308085e-16; + double ans, ax, c, yc, r, t, y, z; + double pk, pkm1, pkm2, qk, qkm1, qkm2; + + /* Compute x**a * exp(-x) / gamma(a) */ + ax = a * log(x) - x - cephes_lgam(a); + if (ax < -MAXLOG) return 0.0; // underflow + ax = exp(ax); + + /* continued fraction */ + y = 1.0 - a; + z = x + y + 1.0; + c = 0.0; + pkm2 = 1.0; + qkm2 = x; + pkm1 = x + 1.0; + qkm1 = z * x; + ans = pkm1/qkm1; + + do { + c += 1.0; + y += 1.0; + z += 2.0; + yc = y * c; + pk = pkm1 * z - pkm2 * yc; + qk = qkm1 * z - qkm2 * yc; + if (qk != 0) { + r = pk/qk; + t = fabs( (ans - r)/r ); + ans = r; + } else { + t = 1.0; + } + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + if (fabs(pk) > BIG) { + pkm2 *= BIGINV; + pkm1 *= BIGINV; + qkm2 *= BIGINV; + qkm1 *= BIGINV; + } + } while( t > MACHEP ); + + return( ans * ax ); +} + + +double cephes_igam(double a, double x) +{ + const double MACHEP = 1.11022302462515654042E-16; // IEEE 2**-53 + const double MAXLOG = 7.09782712893383996843E2; // IEEE log(2**1024) denormalized + double ans, ax, c, r; + + /* Compute x**a * exp(-x) / gamma(a) */ + ax = a * log(x) - x - cephes_lgam(a); + if (ax < -MAXLOG) return 0.0; // underflow + ax = exp(ax); + + /* power series */ + r = a; + c = 1.0; + ans = 1.0; + + do { + r += 1.0; + c *= x/r; + ans += c; + } while (c/ans > MACHEP); + + return ans * ax/a; +} + +/* Logarithm of gamma function */ + + +double cephes_lgam(double x) +{ + const double LOGPI = 1.14472988584940017414; // log(pi) + int sgngam = 1; + double p, q, w, z; + int i; + + if (x < -34.0) { + q = -x; + w = cephes_lgam2(q); + p = floor(q); + if (p == q) return INFINITY; + i = p; + sgngam = ((i&1) == 0) ? -1 : 1; + z = q - p; + if (z > 0.5) { + p += 1.0; + z = p - q; + } + z = q * sin(M_PI * z); + if (z == 0.0) return INFINITY; + z = LOGPI - log(z) - w; + return z; + } else { + return cephes_lgam2(x); + } +} + +double cephes_lgam2(double x) +{ + const double LS2PI = 0.91893853320467274178; // log(sqrt(2*pi)) + const double MAXLGM = 2.556348e305; + int sgngam = 1; + double p, q, u, z; + double A, B, C; + + if (x < 13.0) { + z = 1.0; + p = 0.0; + u = x; + while (u >= 3.0) { + p -= 1.0; + u = x + p; + z *= u; + } + while (u < 2.0) { + if (u == 0.0) return INFINITY; + z /= u; + p += 1.0; + u = x + p; + } + if (z < 0.0) { + sgngam = -1; + z = -z; + } else { + sgngam = 1; + } + if (u == 2.0) return log(z); + p -= 2.0; + x = x + p; + B = (((((( + -1.37825152569120859100E3)*x + -3.88016315134637840924E4)*x + -3.31612992738871184744E5)*x + -1.16237097492762307383E6)*x + -1.72173700820839662146E6)*x + -8.53555664245765465627E5); + C = (((((( + /* 1.00000000000000000000E0)* */x + -3.51815701436523470549E2)*x + -1.70642106651881159223E4)*x + -2.20528590553854454839E5)*x + -1.13933444367982507207E6)*x + -2.53252307177582951285E6)*x + -2.01889141433532773231E6); + p = x * B / C; + return log(z) + p; + } + + if (x > MAXLGM) return sgngam * INFINITY; + + q = (x - 0.5) * log(x) - x + LS2PI; + if (x > 1.0e8) return q; + + p = 1.0/(x*x); + if (x >= 1000.0) { + q += ((7.9365079365079365079365e-4 * p + - 2.7777777777777777777778e-3) * p + + 0.0833333333333333333333) / x; + } else { + A = ((((( + +8.11614167470508450300E-4)*p + -5.95061904284301438324E-4)*p + +7.93650340457716943945E-4)*p + -2.77777777730099687205E-3)*p + +8.33333333333331927722E-2); + q += A / x; + } + return q; +} + +#else // single precision + +// based on cephes/float/igam.c +float gammalnf(float x); +float gammaincf(float a, float x); +float gammainccf(float a, float x); + +float cephes_igamcf(float a, float x); +float cephes_igamf(float a, float x); +float cephes_lgamf(float x); +float cephes_lgam2f(float x); + +// Note: original uses logf, fabsf, floorf and sinf + +float gammaincf(float a, float x) +{ + if ((x <= 0) || ( a <= 0)) return 0.0; + if ((x > 1.0) && (x > a)) return 1.0 - cephes_igamcf(a, x); + return cephes_igamf(a, x); +} + +float gammainccf(float a, float x) +{ + if ((x <= 0) || (a <= 0)) return 1.0; + if ((x < 1.0) || (x < a)) return 1.0 - cephes_igamf(a, x); + return cephes_igamcf(a, x); +} + +float gammalnf(float x) +{ + if (isnan(x)) return(x); + if (!isfinite(x)) return(INFINITY); + return cephes_lgamf(x); +} + +float cephes_igamcf(float a, float x) +{ + const float MAXLOGF = 88.72283905206835; + const float MACHEPF = 5.9604644775390625E-8; + const float BIG = 16777216.; + float ans, c, yc, ax, y, z; + float pk, pkm1, pkm2, qk, qkm1, qkm2; + float r, t; + + ax = a * log(x) - x - cephes_lgamf(a); + if (ax < -MAXLOGF) return 0.0; // underflow + ax = expf(ax); + + /* continued fraction */ + y = 1.0 - a; + z = x + y + 1.0; + c = 0.0; + pkm2 = 1.0; + qkm2 = x; + pkm1 = x + 1.0; + qkm1 = z * x; + ans = pkm1/qkm1; + + do { + c += 1.0; + y += 1.0; + z += 2.0; + yc = y * c; + pk = pkm1 * z - pkm2 * yc; + qk = qkm1 * z - qkm2 * yc; + if (qk != 0) { + r = pk/qk; + t = fabs((ans - r)/r); + ans = r; + } else { + t = 1.0; + } + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + if (fabs(pk) > BIG) { + pkm2 *= MACHEPF; + pkm1 *= MACHEPF; + qkm2 *= MACHEPF; + qkm1 *= MACHEPF; + } + } while (t > MACHEPF); + + return ans * ax; +} + +float cephes_igamf(float a, float x) +{ + const float MAXLOGF = 88.72283905206835; + const float MACHEPF = 5.9604644775390625E-8; + float ans, ax, c, r; + + /* Compute x**a * exp(-x) / gamma(a) */ + ax = a * log(x) - x - cephes_lgamf(a); + if (ax < -MAXLOGF) return 0.0; // underflow + ax = expf(ax); + + /* power series */ + r = a; + c = 1.0; + ans = 1.0; + + do { + r += 1.0; + c *= x/r; + ans += c; + } while (c/ans > MACHEPF); + + return ans * ax/a; +} + +float cephes_lgamf(float x) +{ + const float PIINV = 0.318309886183790671538; + float p, q, w, z; + int i; + int sgngamf; + + if (x < 0.0) { + q = -x; + w = cephes_lgam2f(q); + p = floor(q); + if (p == q) return INFINITY; + i = p; + sgngamf = ((i&1) == 0) ? -1 : 1; + z = q - p; + if (z > 0.5) { + p += 1.0; + z = p - q; + } + z = q * sin(M_PI * z); + if (z == 0.0) return sgngamf * INFINITY; + z = -log(PIINV*z) - w; + return z; + } else { + return cephes_lgam2f(x); + } +} + +float cephes_lgam2f(float x) +{ + const float LS2PI = 0.91893853320467274178; // log(sqrt(2*pi)) + const float MAXLGM = 2.035093e36; + float p, q, z; + float nx, tx; + int direction; + int sgngamf = 1; + + if (x < 6.5) { + direction = 0; + z = 1.0; + tx = x; + nx = 0.0; + if (x >= 1.5) { + while (tx > 2.5) { + nx -= 1.0; + tx = x + nx; + z *=tx; + } + x += nx - 2.0; + } else if (x >= 1.25) { + z *= x; + x -= 1.0; /* x + 1 - 2 */ + direction = 1; + } else if (x >= 0.75) { + x -= 1.0; + /* log gamma(x+1), -.25 < x < .25 */ + // p = x * polevlf( x, C, 7 ); + p = (((((((( + +1.369488127325832E-001)*x + -1.590086327657347E-001)*x + +1.692415923504637E-001)*x + -2.067882815621965E-001)*x + +2.705806208275915E-001)*x + -4.006931650563372E-001)*x + +8.224670749082976E-001)*x + -5.772156501719101E-001)*x; + q = 0.0; + return p + q; + } else { + while (tx < 1.5) { + if (tx == 0.0) return sgngamf * INFINITY; + z *=tx; + nx += 1.0; + tx = x + nx; + } + direction = 1; + x += nx - 2.0; + } + + /* log gamma(x+2), -.5 < x < .5 */ + // p = x * polevlf( x, B, 7 ); + p = (((((((( + +6.055172732649237E-004)*x + -1.311620815545743E-003)*x + +2.863437556468661E-003)*x + -7.366775108654962E-003)*x + +2.058355474821512E-002)*x + -6.735323259371034E-002)*x + +3.224669577325661E-001)*x + +4.227843421859038E-001)*x; + + if (z < 0.0) { + sgngamf = -1; + z = -z; + } else { + sgngamf = 1; + } + q = log(z); + if (direction) q = -q; + return p + q; + } else if (x > MAXLGM) { + return sgngamf * INFINITY; + } else { + /* Note, though an asymptotic formula could be used for x >= 3, + * there is cancellation error in the following if x < 6.5. + */ + q = LS2PI - x; + q += (x - 0.5) * log(x); + + if (x <= 1.0e4) { + z = 1.0/x; + p = z * z; + q += ((6.789774945028216E-004 * p + - 2.769887652139868E-003 ) * p + + 8.333316229807355E-002 ) * z; + } + return q; + } +} + +#endif // !single precision + +#if FLOAT_SIZE>4 +#define sas_gammaln gammaln +#define sas_gammainc gammainc +#define sas_gammaincc gammaincc +#else +#define sas_gammaln gammalnf +#define sas_gammainc gammaincf +#define sas_gammaincc gammainccf +#endif diff -Nru sasmodels-0.98/sasmodels/models/spinodal.py sasmodels-0.99/sasmodels/models/spinodal.py --- sasmodels-0.98/sasmodels/models/spinodal.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/models/spinodal.py 2019-02-10 00:46:47.000000000 +0000 @@ -11,7 +11,20 @@ where $x=q/q_0$, $q_0$ is the peak position, $I_{max}$ is the intensity at $q_0$ (parameterised as the $scale$ parameter), and $B$ is a flat -background. The spinodal wavelength is given by $2\pi/q_0$. +background. The spinodal wavelength, $\Lambda$, is given by $2\pi/q_0$. + +The definition of $I_{max}$ in the literature varies. Hashimoto *et al* (1991) +define it as + +.. math:: + I_{max} = \Lambda^3\Delta\rho^2 + +whereas Meier & Strobl (1987) give + +.. math:: + I_{max} = V_z\Delta\rho^2 + +where $V_z$ is the volume per monomer unit. The exponent $\gamma$ is equal to $d+1$ for off-critical concentration mixtures (smooth interfaces) and $2d$ for critical concentration mixtures @@ -27,14 +40,22 @@ ---------- H. Furukawa. Dynamics-scaling theory for phase-separating unmixing mixtures: -Growth rates of droplets and scaling properties of autocorrelation functions. -Physica A 123,497 (1984). +Growth rates of droplets and scaling properties of autocorrelation functions. +Physica A 123, 497 (1984). + +H. Meier & G. Strobl. Small-Angle X-ray Scattering Study of Spinodal +Decomposition in Polystyrene/Poly(styrene-co-bromostyrene) Blends. +Macromolecules 20, 649-654 (1987). + +T. Hashimoto, M. Takenaka & H. Jinnai. Scattering Studies of Self-Assembling +Processes of Polymer Blends in Spinodal Decomposition. +J. Appl. Cryst. 24, 457-466 (1991). Revision History ---------------- * **Author:** Dirk Honecker **Date:** Oct 7, 2016 -* **Revised:** Steve King **Date:** Sep 7, 2018 +* **Revised:** Steve King **Date:** Oct 25, 2018 """ import numpy as np diff -Nru sasmodels-0.98/sasmodels/model_test.py sasmodels-0.99/sasmodels/model_test.py --- sasmodels-0.98/sasmodels/model_test.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/model_test.py 2019-02-10 00:46:47.000000000 +0000 @@ -46,6 +46,7 @@ import sys import unittest +import traceback try: from StringIO import StringIO @@ -73,7 +74,6 @@ from .kernel import KernelModel # pylint: enable=unused-import - def make_suite(loaders, models): # type: (List[str], List[str]) -> unittest.TestSuite """ @@ -85,7 +85,6 @@ *models* is the list of models to test, or *["all"]* to test all models. """ - ModelTestCase = _hide_model_case_from_nose() suite = unittest.TestSuite() if models[0] in core.KINDS: @@ -94,72 +93,76 @@ else: skip = [] for model_name in models: - if model_name in skip: - continue - model_info = load_model_info(model_name) + if model_name not in skip: + model_info = load_model_info(model_name) + _add_model_to_suite(loaders, suite, model_info) + + return suite - #print('------') - #print('found tests in', model_name) - #print('------') - - # if ispy then use the dll loader to call pykernel - # don't try to call cl kernel since it will not be - # available in some environmentes. - 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 = "%s-python"%model_name - test_method_name = "test_%s_python" % model_info.id +def _add_model_to_suite(loaders, suite, model_info): + ModelTestCase = _hide_model_case_from_nose() + + #print('------') + #print('found tests in', model_name) + #print('------') + + # if ispy then use the dll loader to call pykernel + # don't try to call cl kernel since it will not be + # available in some environmentes. + 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 = "%s-python"%model_info.name + test_method_name = "test_%s_python" % model_info.id + test = ModelTestCase(test_name, model_info, + test_method_name, + platform="dll", # so that + dtype="double", + stash=stash) + suite.addTest(test) + else: # kernel implemented in C + + # test using dll if desired + if 'dll' in loaders or not use_opencl(): + test_name = "%s-dll"%model_info.name + test_method_name = "test_%s_dll" % model_info.id test = ModelTestCase(test_name, model_info, - test_method_name, - platform="dll", # so that - dtype="double", - stash=stash) + test_method_name, + platform="dll", + dtype="double", + stash=stash) suite.addTest(test) - else: # kernel implemented in C - # test using dll if desired - if 'dll' in loaders or not use_opencl(): - test_name = "%s-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 use_opencl(): - test_name = "%s-opencl"%model_name - test_method_name = "test_%s_opencl" % model_info.id - # Using dtype=None so that the models that are only - # correct for double precision are not tested using - # single precision. The choice is determined by the - # presence of *single=False* in the model file. - test = ModelTestCase(test_name, model_info, - test_method_name, - platform="ocl", dtype=None, - stash=stash) - #print("defining", test_name) - suite.addTest(test) + # test using opencl if desired and available + if 'opencl' in loaders and use_opencl(): + test_name = "%s-opencl"%model_info.name + test_method_name = "test_%s_opencl" % model_info.id + # Using dtype=None so that the models that are only + # correct for double precision are not tested using + # single precision. The choice is determined by the + # presence of *single=False* in the model file. + test = ModelTestCase(test_name, model_info, + test_method_name, + platform="ocl", dtype=None, + stash=stash) + #print("defining", test_name) + suite.addTest(test) - return suite def _hide_model_case_from_nose(): # type: () -> type @@ -347,13 +350,25 @@ shift = 10**math.ceil(math.log10(abs(target))) return abs(target-actual)/shift < 1.5*10**-digits -def run_one(model): - # type: (str) -> str +# CRUFT: old interface; should be deprecated and removed +def run_one(model_name): + # msg = "use check_model(model_info) rather than run_one(model_name)" + # warnings.warn(msg, category=DeprecationWarning, stacklevel=2) + try: + model_info = load_model_info(model_name) + except Exception: + output = traceback.format_exc() + return output + + success, output = check_model(model_info) + return output + +def check_model(model_info): + # type: (ModelInfo) -> str """ - Run the tests for a single model, printing the results to stdout. + Run the tests for a single model, capturing the output. - *model* can by a python file, which is handy for checking user defined - plugin models. + Returns success status and the output string. """ # Note that running main() directly did not work from within the # wxPython pycrust console. Instead of the results appearing in the @@ -368,13 +383,8 @@ # Build a test suite containing just the model loaders = ['opencl'] if use_opencl() else ['dll'] - models = [model] - try: - suite = make_suite(loaders, models) - except Exception: - import traceback - stream.writeln(traceback.format_exc()) - return + suite = unittest.TestSuite() + _add_model_to_suite(loaders, suite, model_info) # Warn if there are no user defined tests. # Note: the test suite constructed above only has one test in it, which @@ -389,7 +399,7 @@ # for user tests before running the suite. for test in suite: if not test.info.tests: - stream.writeln("Note: %s has no user defined tests."%model) + stream.writeln("Note: %s has no user defined tests."%model_info.name) break else: stream.writeln("Note: no test suite created --- this should never happen") @@ -405,7 +415,7 @@ output = stream.getvalue() stream.close() - return output + return result.wasSuccessful(), output def main(*models): diff -Nru sasmodels-0.98/sasmodels/multiscat.py sasmodels-0.99/sasmodels/multiscat.py --- sasmodels-0.98/sasmodels/multiscat.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/multiscat.py 2019-02-10 00:46:47.000000000 +0000 @@ -341,12 +341,10 @@ given it will use $n_q = 2^k$ such that $\Delta q < q_\mathrm{min}$. *probability* is related to the expected number of scattering - events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$. As a - hack to allow probability to be a fitted parameter, the "value" - can be a function that takes no parameters and returns the current - value of the probability. *coverage* determines how many scattering - steps to consider. The default is 0.99, which sets $n$ such that - $1 \ldots n$ covers 99% of the Poisson probability mass function. + events in the sample $\lambda$ as $p = 1 - e^{-\lambda}$. + *coverage* determines how many scattering steps to consider. The + default is 0.99, which sets $n$ such that $1 \ldots n$ covers 99% + of the Poisson probability mass function. *is2d* is True then 2D scattering is used, otherwise it accepts and returns 1D scattering. @@ -398,7 +396,7 @@ self.qmax = qmax self.qmin = qmin self.nq = nq - self.probability = probability + self.probability = 0. if probability is None else probability self.coverage = coverage self.is2d = is2d self.window = window @@ -455,6 +453,11 @@ self.Iq = None # type: np.ndarray self.Iqxy = None # type: np.ndarray + # Label probability as a fittable parameter, and give its external name + # Note that the external name must be a valid python identifier, since + # is will be set as an experiment attribute. + self.fittable = {'probability': 'scattering_probability'} + def apply(self, theory): if self.is2d: Iq_calc = theory @@ -462,6 +465,7 @@ Iq_calc = np.interp(self._radius, self.q_calc[0], theory) Iq_calc = Iq_calc.reshape(self.nq, self.nq) + # CRUFT: don't need probability as a function anymore probability = self.probability() if callable(self.probability) else self.probability coverage = self.coverage #t0 = time.time() diff -Nru sasmodels-0.98/sasmodels/sasview_model.py sasmodels-0.99/sasmodels/sasview_model.py --- sasmodels-0.98/sasmodels/sasview_model.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/sasview_model.py 2019-02-10 00:46:47.000000000 +0000 @@ -24,11 +24,13 @@ from . import core from . import custom +from . import kernelcl from . import product from . import generate from . import weights from . import modelinfo from .details import make_kernel_args, dispersion_mesh +from .kernelcl import reset_environment # pylint: disable=unused-import try: @@ -68,6 +70,18 @@ #: has changed since we last reloaded. _CACHED_MODULE = {} # type: Dict[str, "module"] +def reset_environment(): + # type: () -> None + """ + Clear the compute engine context so that the GUI can change devices. + + This removes all compiled kernels, even those that are active on fit + pages, but they will be restored the next time they are needed. + """ + kernelcl.reset_environment() + for model in MODELS.values(): + model._model = None + def find_model(modelname): # type: (str) -> SasviewModelType """ @@ -383,6 +397,12 @@ hidden.add('background') self._model_info.parameters.defaults['background'] = 0. + # Update the parameter lists to exclude any hidden parameters + self.magnetic_params = tuple(pname for pname in self.magnetic_params + if pname not in hidden) + self.orientation_params = tuple(pname for pname in self.orientation_params + if pname not in hidden) + self._persistency_dict = {} self.params = collections.OrderedDict() self.dispersion = collections.OrderedDict() @@ -673,7 +693,11 @@ def _calculate_Iq(self, qx, qy=None): if self._model is None: - self._model = core.build_model(self._model_info) + # Only need one copy of the compiled kernel regardless of how many + # times it is used, so store it in the class. Also, to reset the + # compute engine, need to clear out all existing compiled kernels, + # which is much easier to do if we store them in the class. + self.__class__._model = core.build_model(self._model_info) if qy is not None: q_vectors = [np.asarray(qx), np.asarray(qy)] else: @@ -796,6 +820,16 @@ value = self.params[par.name] return value, [value], [1.0] + @classmethod + def runTests(cls): + """ + Run any tests built into the model and captures the test output. + + Returns success flag and output + """ + from .model_test import check_model + return check_model(cls._model_info) + def test_cylinder(): # type: () -> float """ @@ -883,7 +917,7 @@ def magnetic_demo(): Model = _make_standard_model('sphere') model = Model() - model.setParam('M0:sld', 8) + model.setParam('sld_M0', 8) q = np.linspace(-0.35, 0.35, 500) qx, qy = np.meshgrid(q, q) result = model.calculate_Iq(qx.flatten(), qy.flatten()) diff -Nru sasmodels-0.98/sasmodels/special.py sasmodels-0.99/sasmodels/special.py --- sasmodels-0.98/sasmodels/special.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/sasmodels/special.py 2019-02-10 00:46:47.000000000 +0000 @@ -113,6 +113,18 @@ The standard math function, tgamma(x) is unstable for $x < 1$ on some platforms. + sas_gammaln(x): + log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. + + The standard math function, lgamma(x), is incorrect for single + precision on some platforms. + + sas_gammainc(a, x), sas_gammaincc(a, x): + Incomplete gamma function + sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ + and complementary incomplete gamma function + sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ + sas_erf(x), sas_erfc(x): Error function $\text{sas_erf}(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$ @@ -206,6 +218,9 @@ from numpy import fabs, fmin, fmax, trunc, rint from numpy import pi, nan, inf from scipy.special import gamma as sas_gamma +from scipy.special import gammaln as sas_gammaln +from scipy.special import gammainc as sas_gammainc +from scipy.special import gammaincc as sas_gammaincc from scipy.special import erf as sas_erf from scipy.special import erfc as sas_erfc from scipy.special import j0 as sas_J0 diff -Nru sasmodels-0.98/setup.py sasmodels-0.99/setup.py --- sasmodels-0.98/setup.py 2018-10-03 14:02:40.000000000 +0000 +++ sasmodels-0.99/setup.py 2019-02-10 00:46:47.000000000 +0000 @@ -29,6 +29,11 @@ return version[1:-1] raise RuntimeError("Could not read version from %s/__init__.py"%package) +install_requires = ['numpy', 'scipy'] + +if sys.platform=='win32' or sys.platform=='cygwin': + install_requires.append('tinycc') + setup( name='sasmodels', version=find_version('sasmodels'), @@ -60,12 +65,11 @@ 'sasmodels.models': ['*.c', 'lib/*.c'], 'sasmodels': ['*.c', '*.cl'], }, - install_requires=[ - ], + install_requires=install_requires, extras_require={ + 'full': ['docutils', 'bumps', 'matplotlib'], + 'server': ['bumps'], 'OpenCL': ["pyopencl"], - 'Bumps': ["bumps"], - 'TinyCC': ["tinycc"], }, build_requires=['setuptools'], test_requires=['pytest'],