diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index af31633..d8b8f4c 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -10,4 +10,3 @@ Release notes
**Added:**
* Add source code from ``diffpy.srfit.pdf`` to ``diffpy.cmipdf``.
-
diff --git a/CODE-OF-CONDUCT.rst b/CODE-OF-CONDUCT.rst
index e8199ca..25fafe2 100644
--- a/CODE-OF-CONDUCT.rst
+++ b/CODE-OF-CONDUCT.rst
@@ -67,7 +67,7 @@ Enforcement
Instances of abusive, harassing, or otherwise unacceptable behavior may be
reported to the community leaders responsible for enforcement at
-sb2896@columbia.edu. All complaints will be reviewed and investigated promptly and fairly.
+sbillinge@ucsb.edu. All complaints will be reviewed and investigated promptly and fairly.
All community leaders are obligated to respect the privacy and security of the
reporter of any incident.
diff --git a/LICENSE_DANSE.rst b/LICENSE_DANSE.rst
new file mode 100644
index 0000000..9e65c74
--- /dev/null
+++ b/LICENSE_DANSE.rst
@@ -0,0 +1,34 @@
+This program is part of the DiffPy and DANSE open-source projects at Columbia
+University and is available subject to the conditions and terms laid out below.
+
+Copyright (c) 2008-2011, The Trustees of Columbia University in
+the City of New York. All rights reserved.
+
+For more information please visit the diffpy web-page at
+ http://www.diffpy.org
+or email Prof. Simon Billinge at sbillinge@ucsb.edu.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright notice, this
+ list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above copyright notice,
+ this list of conditions and the following disclaimer in the documentation
+ and/or other materials provided with the distribution.
+
+ * Neither the names of COLUMBIA UNIVERSITY, MICHIGAN STATE UNIVERSITY nor the
+ names of their contributors may be used to endorse or promote products
+ derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.rst b/README.rst
index 44de95f..95051b3 100644
--- a/README.rst
+++ b/README.rst
@@ -129,7 +129,7 @@ Before contributing, please read our `Code of Conduct `_ or email Simon Billinge at sb2896@columbia.edu.
+For more information on diffpy.cmipdf please visit the project `web-page `_ or email Simon Billinge at sbillinge@ucsb.edu.
Acknowledgements
----------------
diff --git a/cookiecutter.json b/cookiecutter.json
index 95bf05e..f5a6d11 100644
--- a/cookiecutter.json
+++ b/cookiecutter.json
@@ -1,6 +1,6 @@
{
"maintainer_name": "Simon Billinge",
- "maintainer_email": "sb2896@columbia.edu",
+ "maintainer_email": "sbillinge@ucsb.edu",
"maintainer_github_username": "sbillinge",
"contributors": "Caden Myers, Simon Billinge, and members of the Billinge group",
"license_holders": "Simon Billinge",
diff --git a/news/fix-cfs.rst b/news/fix-cfs.rst
new file mode 100644
index 0000000..4c23e04
--- /dev/null
+++ b/news/fix-cfs.rst
@@ -0,0 +1,24 @@
+**Added:**
+
+*
+
+**Changed:**
+
+* Rename the parameters of the characteristic functions ``spherical_particle``, ``spheroidal_particle``, ``lognormal_spherical_distribution``, ``sheet_particle``, and ``spherical_shell`` to match ``diffpy.srfit.pdf.characteristicfunctions``.
+* Issue a ``RuntimeWarning`` once per process when a characteristic function is given non-physical input and returns zero.
+
+**Deprecated:**
+
+*
+
+**Removed:**
+
+* Remove the non-functional ``SASCF`` calculator from ``characteristicfunctions``. SAS characteristic functions now live in ``diffpy.cmisas``.
+
+**Fixed:**
+
+*
+
+**Security:**
+
+*
diff --git a/pyproject.toml b/pyproject.toml
index a1024d2..9419418 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -6,10 +6,10 @@ build-backend = "setuptools.build_meta"
name = "diffpy.cmipdf"
dynamic=['version', 'dependencies']
authors = [
- { name="Simon Billinge", email="sb2896@columbia.edu" },
+ { name="Simon Billinge", email="sbillinge@ucsb.edu" },
]
maintainers = [
- { name="Simon Billinge", email="sb2896@columbia.edu" },
+ { name="Simon Billinge", email="sbillinge@ucsb.edu" },
]
description = "The code of the PDF pack of the diffpy.cmi package."
keywords = ['diffraction', 'PDF', 'X-ray', 'neutron']
diff --git a/src/diffpy/cmipdf/characteristicfunctions.py b/src/diffpy/cmipdf/characteristicfunctions.py
index bebf880..90eab31 100644
--- a/src/diffpy/cmipdf/characteristicfunctions.py
+++ b/src/diffpy/cmipdf/characteristicfunctions.py
@@ -31,104 +31,143 @@
"lognormal_spherical_distribution",
"sheet_particle",
"spherical_shell",
- "SASCF",
]
+import warnings
+
import numpy
from numpy import arctan as atan
from numpy import arctanh as atanh
-from numpy import ceil, exp, log, log2, pi, sign, sqrt
-from numpy.fft import fftfreq, ifft
+from numpy import exp, log, sign, sqrt
from scipy.special import erf
-from diffpy.srfit.fitbase.calculator import Calculator
-
-
-def spherical_particle(radial_dist, p_diameter):
- """Spherical nanoparticle characteristic function.
+# The record of non-physical input warnings already issued in this process,
+# keyed by (function name, requirement). See _warn_non_physical.
+_warned_non_physical = set()
- Parameters
- ----------
- radial_dist : float or array-like
- The distance of interaction.
- p_diameter : float
- The particle diameter.
+def _warn_non_physical(function_name, requirement):
+ """Warn that non-physical input produced a zero characteristic
+ function.
- From Kodama et al., Acta Cryst. A, 62, 444-453
- (converted from radius to diameter)
+ A refinement evaluates a characteristic function on every iteration,
+ so the warning is issued only once per process for each distinct
+ requirement. The message must not include the offending value, since
+ the `warnings` registry keys on the message text and a value that
+ drifts from iteration to iteration would defeat the suppression.
"""
- f = numpy.zeros(numpy.shape(radial_dist), dtype=float)
- if p_diameter > 0:
- x = numpy.array(radial_dist, dtype=float) / p_diameter
- inside = x < 1.0
- xin = x[inside]
- f[inside] = 1.0 - 1.5 * xin + 0.5 * xin * xin * xin
- return f
+ key = (function_name, requirement)
+ if key in _warned_non_physical:
+ return
+ _warned_non_physical.add(key)
+ warnings.warn(
+ f"In '{function_name}', {requirement}. The characteristic function "
+ "was set to zero for every r, which flattens the fit residual so a "
+ "refinement cannot recover on its own. Please use a physically "
+ "meaningful starting value, or keep the parameter in range with "
+ "FitRecipe.add_soft_bounds.",
+ RuntimeWarning,
+ stacklevel=3,
+ )
-def spheroidal_particle(radial_dist, r_equatorial, r_polar):
- """Spheroidal characteristic function specified using radii.
+def spherical_particle(r, particle_diameter):
+ """Compute the spherical nanoparticle characteristic function.
- Spheroid with radii (r_equatorial, r_equatorial, r_polar)
+ From Kodama et al., Acta Cryst. A, 62, 444-453 (converted from
+ radius to diameter).
Parameters
----------
- radial_dist : float or array-like
+ r : array_like
The distance of interaction.
- r_polar : float
- The polar radius of the spheroid.
- r_equatorial
- The equatorial radius of the spheroid.
-
+ particle_diameter : float
+ The particle diameter.
- Note
- ----
- - `r_equatorial < r_polar` equates to a prolate spheroid
- - `r_equatorial > r_polar` equates to a oblate spheroid
- - `r_equatorial == r_polar` is a sphere
+ Returns
+ -------
+ numpy.ndarray
+ The characteristic function values evaluated at `r`.
+
+ Warns
+ -----
+ RuntimeWarning
+ If `particle_diameter` is not positive, in which case the
+ characteristic function is zero everywhere. The warning is issued
+ only once per process.
"""
- d_equatorial = 2.0 * r_equatorial
- pelpt = 1.0 * r_polar / r_equatorial
- return _calculate_spheroidal_cf(radial_dist, d_equatorial, pelpt)
+ characteristic_function = numpy.zeros(numpy.shape(r), dtype=float)
+ if particle_diameter <= 0:
+ _warn_non_physical(
+ "spherical_particle", "'particle_diameter' must be positive"
+ )
+ return characteristic_function
+ scaled_r = numpy.array(r, dtype=float) / particle_diameter
+ inside = scaled_r < 1.0
+ scaled_r_inside = scaled_r[inside]
+ characteristic_function[inside] = (
+ 1.0 - 1.5 * scaled_r_inside + 0.5 * scaled_r_inside**3
+ )
+ return characteristic_function
-def _calculate_spheroidal_cf(r, d_equatorial, axis_ratio):
- """Calculate the spheroidal nanoparticle characteristic function.
+def spheroidal_particle(r, equatorial_radius, polar_radius):
+ """Compute the spheroidal nanoparticle characteristic function.
- Form factor for ellipsoid with radii
- (d_equatorial/2, d_equatorial/2, axis_ratio*d_equatorial/2)
+ Spheroid with radii (equatorial_radius, equatorial_radius,
+ polar_radius). ``equatorial_radius < polar_radius`` equates to a
+ prolate spheroid, ``equatorial_radius > polar_radius`` equates to
+ an oblate spheroid, and ``equatorial_radius == polar_radius`` is a
+ sphere. From Lei et al., Phys. Rev. B, 80, 024118 (2009).
Parameters
----------
- r : float or array-like
- The distance of interaction
- d_equatorial : float
- The equatorial diameter
- axis_ratio : float
- The ratio of axis lengths
-
-
- From Lei et al., Phys. Rev. B, 80, 024118 (2009)
+ r : array_like
+ The distance of interaction.
+ equatorial_radius : float
+ The equatorial radius.
+ polar_radius : float
+ The polar radius.
+
+ Returns
+ -------
+ numpy.ndarray
+ The characteristic function values evaluated at `r`.
+
+ Warns
+ -----
+ RuntimeWarning
+ If `equatorial_radius` or `polar_radius` is not positive, in which
+ case the characteristic function is zero everywhere. The warning is
+ issued only once per process.
"""
- pelpt = 1.0 * axis_ratio
+ if equatorial_radius <= 0:
+ _warn_non_physical(
+ "spheroidal_particle", "'equatorial_radius' must be positive"
+ )
+ return numpy.zeros(numpy.shape(r), dtype=float)
+ if polar_radius <= 0:
+ _warn_non_physical(
+ "spheroidal_particle", "'polar_radius' must be positive"
+ )
+ return numpy.zeros(numpy.shape(r), dtype=float)
- if d_equatorial <= 0 or pelpt <= 0:
- return numpy.zeros_like(r)
+ particle_diameter = 2.0 * equatorial_radius
+ axis_ratio = 1.0 * polar_radius / equatorial_radius
# to simplify the equations
- v = pelpt
- d = 1.0 * d_equatorial
+ v = axis_ratio
+ d = particle_diameter
d2 = d * d
v2 = v * v
if v == 1:
- return spherical_particle(r, d_equatorial)
+ return spherical_particle(r, particle_diameter)
rx = r
if v < 1:
- r = rx[rx <= v * d_equatorial]
+ r = rx[rx <= v * d]
r2 = r * r
f1 = (
1
@@ -142,7 +181,7 @@ def _calculate_spheroidal_cf(r, d_equatorial, axis_ratio):
* atanh(sqrt(1 - v2))
)
- r = rx[numpy.logical_and(rx > v * d_equatorial, rx <= d_equatorial)]
+ r = rx[numpy.logical_and(rx > v * d, rx <= d)]
r2 = r * r
f2 = (
(
@@ -157,14 +196,14 @@ def _calculate_spheroidal_cf(r, d_equatorial, axis_ratio):
/ sqrt(1 - v2)
)
- r = rx[rx > d_equatorial]
+ r = rx[rx > d]
f3 = numpy.zeros_like(r)
f = numpy.concatenate((f1, f2, f3))
elif v > 1:
- r = rx[rx <= d_equatorial]
+ r = rx[rx <= d]
r2 = r * r
f1 = (
1
@@ -178,7 +217,7 @@ def _calculate_spheroidal_cf(r, d_equatorial, axis_ratio):
* atan(sqrt(v2 - 1))
)
- r = rx[numpy.logical_and(rx > d_equatorial, rx <= v * d_equatorial)]
+ r = rx[numpy.logical_and(rx > d, rx <= v * d)]
r2 = r * r
f2 = (
1
@@ -198,7 +237,7 @@ def _calculate_spheroidal_cf(r, d_equatorial, axis_ratio):
* (atan(sqrt(v2 - 1)) - atan(sqrt(r2 / d2 - 1)))
)
- r = rx[rx > v * d_equatorial]
+ r = rx[rx > v * d]
f3 = numpy.zeros_like(r)
f = numpy.concatenate((f1, f2, f3))
@@ -206,26 +245,16 @@ def _calculate_spheroidal_cf(r, d_equatorial, axis_ratio):
return f
-def lognormal_spherical_distribution(radial_dist, p_diameter, p_sigma):
- """Spherical nanoparticle characteristic function with lognormal
- size distribution.
+def lognormal_spherical_distribution(
+ r, particle_diameter, particle_diameter_sigma
+):
+ """Compute the spherical nanoparticle characteristic function with
+ lognormal size distribution.
- Parameters
- ----------
- radial_dist : float or array-like
- The distance of interaction.
- p_diameter : float
- The mean particle diameter.
- p_sigma : float
- The log-normal width of the particle diameter.
-
-
- Here, radial_dist is the independent variable, mu is the mean of the
- distribution
- (not of the particle size), and s is the width of the distribution.
- This is
- the characteristic function for the lognormal distribution of particle
- diameter::
+ Here, r is the independent variable, mu is the mean of the
+ distribution (not of the particle size), and s is the width of
+ the distribution. This is the characteristic function for the
+ lognormal distribution of particle diameter::
F(r, mu, s) = (
0.5*Erfc((-mu-3*s^2+Log(r))/(sqrt(2)*s))
@@ -234,109 +263,174 @@ def lognormal_spherical_distribution(radial_dist, p_diameter, p_sigma):
)
The expectation value of the distribution gives the average particle
- diameter, p_diameter. The variance of the distribution gives p_sigma^2.
- mu and s can be expressed in terms of these as::
+ diameter, particle_diameter. The variance of the distribution gives
+ particle_diameter_sigma^2. mu and s can be expressed in terms of these
+ as::
- s^2 = log((p_sigma/p_diameter)^2 + 1)
- mu = log(p_diameter) - s^2/2
+ s^2 = log((particle_diameter_sigma/particle_diameter)^2 + 1)
+ mu = log(particle_diameter) - s^2/2
- Source unknown
+ Source unknown.
+
+ Parameters
+ ----------
+ r : array_like
+ The distance of interaction.
+ particle_diameter : float
+ The mean particle diameter.
+ particle_diameter_sigma : float
+ The log-normal width of the particle diameter.
+
+ Returns
+ -------
+ numpy.ndarray
+ The characteristic function values evaluated at `r`.
+
+ Warns
+ -----
+ RuntimeWarning
+ If `particle_diameter` is not positive, if `particle_diameter_sigma`
+ is negative, or if `particle_diameter` is too small for
+ `particle_diameter_sigma` for the closed form to hold, in which case
+ the characteristic function is zero everywhere. The warning is issued
+ only once per process. A `particle_diameter_sigma` of zero is the
+ sphere limit rather than an error and does not warn.
"""
- if p_diameter <= 0:
- return numpy.zeros_like(radial_dist)
- if p_sigma <= 0:
- return spherical_particle(radial_dist, p_diameter)
+ if particle_diameter <= 0:
+ _warn_non_physical(
+ "lognormal_spherical_distribution",
+ "'particle_diameter' must be positive",
+ )
+ return numpy.zeros(numpy.shape(r), dtype=float)
+ if particle_diameter_sigma < 0:
+ _warn_non_physical(
+ "lognormal_spherical_distribution",
+ "'particle_diameter_sigma' must not be negative",
+ )
+ return numpy.zeros(numpy.shape(r), dtype=float)
+ if particle_diameter_sigma == 0:
+ return spherical_particle(r, particle_diameter)
sqrt2 = sqrt(2.0)
- s = sqrt(log(p_sigma * p_sigma / (1.0 * p_diameter * p_diameter) + 1))
- mu = log(p_diameter) - s * s / 2
+ s = sqrt(
+ log(
+ particle_diameter_sigma
+ * particle_diameter_sigma
+ / (1.0 * particle_diameter * particle_diameter)
+ + 1
+ )
+ )
+ mu = log(particle_diameter) - s * s / 2
if mu < 0:
- return numpy.zeros_like(radial_dist)
+ _warn_non_physical(
+ "lognormal_spherical_distribution",
+ "'particle_diameter' is too small for 'particle_diameter_sigma'",
+ )
+ return numpy.zeros(numpy.shape(r), dtype=float)
return (
- 0.5 * erfc((-mu - 3 * s * s + log(radial_dist)) / (sqrt2 * s))
+ 0.5 * erfc((-mu - 3 * s * s + log(r)) / (sqrt2 * s))
+ 0.25
- * radial_dist
- * radial_dist
- * radial_dist
- * erfc((-mu + log(radial_dist)) / (sqrt2 * s))
+ * r
+ * r
+ * r
+ * erfc((-mu + log(r)) / (sqrt2 * s))
* exp(-3 * mu - 4.5 * s * s)
- 0.75
- * radial_dist
- * erfc((-mu - 2 * s * s + log(radial_dist)) / (sqrt2 * s))
+ * r
+ * erfc((-mu - 2 * s * s + log(r)) / (sqrt2 * s))
* exp(-mu - 2.5 * s * s)
)
-def sheet_particle(r, thickness):
- """Nanosheet characteristic function.
+def sheet_particle(r, sheet_thickness):
+ """Compute the nanosheet characteristic function.
+
+ From Kodama et al., Acta Cryst. A, 62, 444-453.
Parameters
----------
- r: float or array-like
+ r : array_like
The distance of interaction.
- thickness : float
- The thickness of nanosheet.
-
-
- From Kodama et al., Acta Cryst. A, 62, 444-453
+ sheet_thickness : float
+ The thickness of the nanosheet.
+
+ Returns
+ -------
+ numpy.ndarray
+ The characteristic function values evaluated at `r`.
+
+ Warns
+ -----
+ RuntimeWarning
+ If `sheet_thickness` is not positive, in which case the
+ characteristic function is zero everywhere. The warning is issued
+ only once per process.
"""
- # handle zero or negative thickness. make it work for scalars and arrays.
- if thickness <= 0:
- return 0 * thickness
+ if sheet_thickness <= 0:
+ _warn_non_physical(
+ "sheet_particle", "'sheet_thickness' must be positive"
+ )
+ if numpy.isscalar(r):
+ return 0.0
+ return numpy.zeros(numpy.shape(r), dtype=float)
# process scalar r
if numpy.isscalar(r):
- rv = 1 - 0.5 * r / thickness if r < thickness else 0.5 * thickness / r
+ rv = (
+ 1 - 0.5 * r / sheet_thickness
+ if r < sheet_thickness
+ else 0.5 * sheet_thickness / r
+ )
return rv
# handle array-type r
- ra = numpy.asarray(r)
- lo = ra < thickness
- hi = ~lo
- f = numpy.empty_like(ra, dtype=float)
- f[lo] = 1 - 0.5 * ra[lo] / thickness
- f[hi] = 0.5 * thickness / ra[hi]
- return f
+ r_array = numpy.asarray(r)
+ inside = r_array < sheet_thickness
+ outside = ~inside
+ characteristic_function = numpy.empty_like(r_array, dtype=float)
+ characteristic_function[inside] = (
+ 1 - 0.5 * r_array[inside] / sheet_thickness
+ )
+ characteristic_function[outside] = 0.5 * sheet_thickness / r_array[outside]
+ return characteristic_function
+
+def spherical_shell(r, radius, thickness):
+ """Compute the spherical shell characteristic function.
-def spherical_shell(radial_dist, inner_radius, thickness):
- """Spherical shell characteristic function.
+ The outer radius equals ``radius + thickness``. From Lei et al.,
+ Phys. Rev. B, 80, 024118 (2009).
Parameters
----------
- radial_dist : float or array-like
+ r : array_like
The distance of interaction.
- inner_radius : float
- The inner radius of the shell.
+ radius : float
+ The inner radius.
thickness : float
- The thickness of shell.
-
-
- outer radius = inner_radius + thickness
-
- From Lei et al., Phys. Rev. B, 80, 024118 (2009)
+ The thickness of the shell.
+
+ Returns
+ -------
+ numpy.ndarray
+ The characteristic function values evaluated at `r`.
+
+ Warns
+ -----
+ RuntimeWarning
+ If `thickness` is not positive or `radius` is negative, in which case
+ the characteristic function is zero everywhere. The warning is issued
+ only once per process. A `radius` of zero is the solid sphere limit
+ rather than an error and does not warn.
"""
- d = 1.0 * thickness
- a = 1.0 * inner_radius + d / 2.0
- return _calculate_shell_cf(radial_dist, a, d)
-
-
-def _calculate_shell_cf(r, a, delta):
- """Spherical shell characteristic function.
-
- Parameters
- ----------
- a
- Central radius
- delta
- Thickness of shell
-
- outer radius = a + thickness/2
-
+ if thickness <= 0:
+ _warn_non_physical("spherical_shell", "'thickness' must be positive")
+ return numpy.zeros(numpy.shape(r), dtype=float)
+ if radius < 0:
+ _warn_non_physical("spherical_shell", "'radius' must not be negative")
+ return numpy.zeros(numpy.shape(r), dtype=float)
- From Lei et al., Phys. Rev. B, 80, 024118 (2009)
- """
- a = 1.0 * a
- d = 1.0 * delta
+ d = 1.0 * thickness
+ a = 1.0 * radius + d / 2.0
a2 = a**2
d2 = d**2
dmr = d - r
@@ -365,126 +459,19 @@ def _calculate_shell_cf(r, a, delta):
return f
-class SASCF(Calculator):
- """Calculator class for characteristic functions from sas-models.
-
- This class wraps a sas.models.BaseModel to calculate I(Q) related to
- nanoparticle shape. This I(Q) is inverted to f(r) according to:
- f(r) = 1 / (4 pi r) * SINFT(I(Q)),
- where "SINFT" represents the sine Fourier transform.
+def erfc(x):
+ """Compute the complementary error function.
- Attributes
+ Parameters
----------
- _model
- BaseModel object this adapts.
- Managed Parameters
- These depend on the parameters of the BaseModel object held by _model.
- They are created from the 'params' attribute of the BaseModel. If a
- dispersion is set for the BaseModel, the dispersion "width" will be
- accessible under "_width", where is the name a
- parameter adjusted by dispersion.
- """
-
- def __init__(self, name, model):
- """Initialize the generator.
+ x : array_like
+ The input value.
- Parameters
- ----------
- name
- A name for the SASCF
- model
- SASModel object this adapts.
- """
- Calculator.__init__(self, name)
-
- self._model = model
-
- from diffpy.srfit.sas.sasparameter import SASParameter
-
- # Wrap normal parameters
- for parname in model.params:
- par = SASParameter(parname, model)
- self.add_parameter(par)
-
- # Wrap dispersion parameters
- for parname in model.dispersion:
- name = parname + "_width"
- parname += ".width"
- par = SASParameter(name, model, parname)
- self.add_parameter(par)
-
- return
-
- def __call__(self, r):
- """Calculate the characteristic function from the transform of
- the BaseModel."""
-
- # Determine q-values.
- # We want very fine r-spacing so we can properly normalize f(r). This
- # equates to having a large qmax so that the Fourier transform is
- # finely spaced. We also want the calculation to be fast, so we pick
- # qmax such that the number of q-points is a power of 2. This allows us
- # to use the fft.
- #
- # The initial dr is somewhat arbitrary, but using dr = 0.01 allows for
- # the f(r) calculated from a particle of diameter 50, over r =
- # arange(1, 60, 0.1) to agree with the spherical_particle with
- # Rw < 1e-4%.
- #
- # We also have to make a q-spacing small enough to compute out to at
- # least the size of the signal.
- raise NotImplementedError(
- "As of release 3.2.0, SAS characteristic functions are not working"
- + " but we hope to have them working again in a future release."
- )
-
- dr = min(0.01, r[1] - r[0])
- ed = 2 * self._model.calculate_ER()
-
- # Check for nans. If we find any, then return zeros.
- if numpy.isnan(ed).any():
- y = numpy.zeros_like(r)
- return y
-
- rmax = max(ed, 2 * r[-1])
- dq = pi / rmax
- qmax = pi / dr
- numpoints = int(2 ** (ceil(log2(qmax / dq))))
- qmax = dq * numpoints
-
- # Calculate F(q) = q * I(q) from model
- q = fftfreq(int(qmax / dq)) * qmax
- fq = q * self._model.evalDistribution(q)
-
- # Calculate g(r) and the effective r-points
- rp = fftfreq(numpoints) * 2 * pi / dq
- # Note sine transform = imaginary part of ifft
- gr = ifft(fq).imag
-
- # Calculate full-fr for normalization
- assert rp[0] == 0.0
- frp = numpy.zeros_like(gr)
- frp[1:] = gr[1:] / rp[1:]
-
- # Inerpolate onto requested grid, do not use data after jump in rp
- assert numpoints % 2 == 0
- nhalf = numpoints / 2
- fr = numpy.interp(r, rp[:nhalf], gr[:nhalf])
- vmask = r != 0
- fr[vmask] /= r[vmask]
-
- # Normalize. We approximate fr[0] by using the fact that f(r) is linear
- # at low r. By definition, fr[0] should equal 1.
- fr0 = 2 * frp[2] - frp[1]
- fr /= fr0
-
- # Fix potential divide-by-zero issue, fr is 1 at r == 0
- fr[~vmask] = 1
-
- return fr
-
-
-def erfc(x):
+ Returns
+ -------
+ numpy.ndarray
+ The complementary error function evaluated at `x`.
+ """
return 1.0 - erf(x)
diff --git a/tests/conftest.py b/tests/conftest.py
index e4544a1..0227883 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -7,6 +7,22 @@
from diffpy.structure import PDFFitStructure
+@pytest.fixture
+def reset_characteristic_function_warnings():
+ """Clear the record of already-issued non-physical input warnings.
+
+ The characteristic functions warn only once per process for each
+ distinct non-physical input so that a refinement loop does not flood
+ the user with repeated messages. Tests that assert on the warning
+ must start and finish with a clean record.
+ """
+ from diffpy.cmipdf.characteristicfunctions import _warned_non_physical
+
+ _warned_non_physical.clear()
+ yield
+ _warned_non_physical.clear()
+
+
@pytest.fixture
def user_filesystem(tmp_path):
base_dir = Path(tmp_path)
diff --git a/tests/test_characteristicfunctions.py b/tests/test_characteristicfunctions.py
new file mode 100644
index 0000000..c194194
--- /dev/null
+++ b/tests/test_characteristicfunctions.py
@@ -0,0 +1,416 @@
+"""Tests for the nanoparticle characteristic functions."""
+
+# FIXME: remove this line when `docformatter` fixes the blank line bug
+
+
+import re
+import unittest
+import warnings
+
+import numpy
+import numpy.testing as npt
+import pytest
+
+import diffpy.cmipdf.characteristicfunctions as cf
+
+# ----------------------------------------------------------------------------
+# The characteristic functions below compute the attenuation of the PDF due
+# to a finite nanoparticle shape. Each has a known value of 1 at r = 0 (the
+# particle always overlaps itself at zero distance) and drops to 0 once r
+# reaches the largest distance spanned by the shape.
+
+
+@pytest.mark.parametrize(
+ "input_r, input_particle_diameter, expected_characteristic_function",
+ [
+ # C1: r ranges over values inside and outside the particle.
+ # Expected: the polynomial form 1 - 1.5x + 0.5x^3 (x = r /
+ # particle_diameter) applies inside the particle, and the function
+ # is zero at and beyond the particle diameter.
+ (
+ numpy.array([0.0, 5.0, 10.0, 15.0]),
+ 10.0,
+ [1.0, 0.3125, 0.0, 0.0],
+ ),
+ ],
+)
+def test_spherical_particle(
+ input_r, input_particle_diameter, expected_characteristic_function
+):
+ actual_characteristic_function = cf.spherical_particle(
+ input_r, input_particle_diameter
+ )
+ npt.assert_allclose(
+ actual_characteristic_function, expected_characteristic_function
+ )
+
+
+@pytest.mark.parametrize(
+ "input_r, input_equatorial_radius, input_polar_radius, "
+ "expected_characteristic_function",
+ [
+ # C1: equatorial_radius equals polar_radius, i.e. a sphere.
+ # Expected: matches spherical_particle for the equivalent diameter.
+ (
+ numpy.linspace(0, 25, 6),
+ 10.0,
+ 10.0,
+ cf.spherical_particle(numpy.linspace(0, 25, 6), 20.0),
+ ),
+ # C2: polar_radius > equatorial_radius, i.e. a prolate spheroid.
+ # Expected: values from the Lei et al. characteristic function.
+ (
+ numpy.array([0.0, 10.0, 20.0, 30.0, 40.0, 50.0]),
+ 10.0,
+ 20.0,
+ [
+ 1.0,
+ 0.40106264900760513,
+ 0.054200238412168256,
+ 0.003650768214115807,
+ 0.0,
+ 0.0,
+ ],
+ ),
+ # C3: polar_radius < equatorial_radius, i.e. an oblate spheroid.
+ # Expected: values from the Lei et al. characteristic function.
+ (
+ numpy.array([0.0, 5.0, 10.0, 15.0, 20.0, 25.0]),
+ 20.0,
+ 10.0,
+ [
+ 1.0,
+ 0.744181556741916,
+ 0.5061470768546105,
+ 0.30368052370886184,
+ 0.15456586067544859,
+ 0.0675583474738014,
+ ],
+ ),
+ ],
+)
+def test_spheroidal_particle(
+ input_r,
+ input_equatorial_radius,
+ input_polar_radius,
+ expected_characteristic_function,
+):
+ actual_characteristic_function = cf.spheroidal_particle(
+ input_r, input_equatorial_radius, input_polar_radius
+ )
+ npt.assert_allclose(
+ actual_characteristic_function, expected_characteristic_function
+ )
+
+
+@pytest.mark.parametrize(
+ "input_r, input_particle_diameter, input_particle_diameter_sigma, "
+ "expected_characteristic_function",
+ [
+ # C1: the lognormal distribution has no width.
+ # Expected: matches spherical_particle for the same diameter.
+ (
+ numpy.array([0.0, 5.0, 10.0]),
+ 10.0,
+ 0.0,
+ cf.spherical_particle(numpy.array([0.0, 5.0, 10.0]), 10.0),
+ ),
+ # C2: a genuine lognormal size distribution is applied.
+ # Expected: values from the closed-form lognormal-averaged
+ # spherical characteristic function.
+ (
+ numpy.array([1.0, 5.0, 10.0, 15.0]),
+ 10.0,
+ 2.0,
+ [
+ 0.8617610662266728,
+ 0.362144891493816,
+ 0.039068792680198694,
+ 0.0009056141323687261,
+ ],
+ ),
+ ],
+)
+def test_lognormal_spherical_distribution(
+ input_r,
+ input_particle_diameter,
+ input_particle_diameter_sigma,
+ expected_characteristic_function,
+):
+ actual_characteristic_function = cf.lognormal_spherical_distribution(
+ input_r, input_particle_diameter, input_particle_diameter_sigma
+ )
+ npt.assert_allclose(
+ actual_characteristic_function, expected_characteristic_function
+ )
+
+
+@pytest.mark.parametrize(
+ "input_r, input_sheet_thickness, expected_characteristic_function",
+ [
+ # C1: r is an array with values inside and outside the sheet
+ # thickness.
+ # Expected: 1 - 0.5*r/thickness inside the sheet, and
+ # 0.5*thickness/r beyond it.
+ (
+ numpy.array([0.0, 2.0, 4.0, 8.0]),
+ 4.0,
+ [1.0, 0.75, 0.5, 0.25],
+ ),
+ # C2: r is a scalar inside the sheet thickness.
+ # Expected: 1 - 0.5*r/thickness.
+ (2.0, 4.0, 0.75),
+ ],
+)
+def test_sheet_particle(
+ input_r, input_sheet_thickness, expected_characteristic_function
+):
+ actual_characteristic_function = cf.sheet_particle(
+ input_r, input_sheet_thickness
+ )
+ npt.assert_allclose(
+ actual_characteristic_function, expected_characteristic_function
+ )
+
+
+@pytest.mark.parametrize(
+ "input_r, input_radius, input_thickness, "
+ "expected_characteristic_function",
+ [
+ # C1: r ranges from the shell center out past the outer radius.
+ # Expected: values from the Lei et al. shell characteristic
+ # function; the function is 1 at r = 0 and 0 once r exceeds the
+ # outer radius (radius + thickness).
+ (
+ numpy.array([0.0, 5.0, 12.5, 15.0, 20.0, 30.0]),
+ 10.0,
+ 5.0,
+ [
+ 1.0,
+ 0.4934210526315789,
+ 0.19736842105263158,
+ 0.16447368421052633,
+ 0.12335526315789473,
+ 0.0,
+ ],
+ ),
+ ],
+)
+def test_spherical_shell(
+ input_r, input_radius, input_thickness, expected_characteristic_function
+):
+ actual_characteristic_function = cf.spherical_shell(
+ input_r, input_radius, input_thickness
+ )
+ npt.assert_allclose(
+ actual_characteristic_function, expected_characteristic_function
+ )
+
+
+# ----------------------------------------------------------------------------
+# A refinement can step a shape parameter into a region that has no physical
+# meaning, such as a negative diameter or thickness. The characteristic
+# functions must not raise there, because that would abort a refinement that
+# would otherwise recover on its own, so they return zero for every r instead.
+# A zero return is flat in the shape parameter and so stalls the optimizer,
+# which the user cannot see from the fit output alone; each function therefore
+# also warns. The warning is issued only once per process for each distinct
+# problem so that a refinement loop does not bury the user in repeated
+# messages.
+#
+# Every such warning states which requirement was violated and then closes
+# with the same explanation of the consequence and of how to avoid it. Only
+# the requirement varies between cases.
+
+
+@pytest.mark.parametrize(
+ "function_name, input_args, expected_requirement",
+ [
+ # C1: the sphere diameter is zero.
+ # Expected: zero everywhere, warning that the diameter must be
+ # positive.
+ (
+ "spherical_particle",
+ (numpy.array([0.0, 5.0]), 0.0),
+ "'particle_diameter' must be positive",
+ ),
+ # C2: the sphere diameter is negative.
+ # Expected: zero everywhere, warning that the diameter must be
+ # positive.
+ (
+ "spherical_particle",
+ (numpy.array([1.0, 5.0, 10.0]), -10.0),
+ "'particle_diameter' must be positive",
+ ),
+ # C3: the spheroid equatorial radius is zero, which the axis ratio
+ # would otherwise be divided by.
+ # Expected: zero everywhere, warning that the equatorial radius must
+ # be positive, and no ZeroDivisionError.
+ (
+ "spheroidal_particle",
+ (numpy.array([1.0, 5.0]), 0.0, 5.0),
+ "'equatorial_radius' must be positive",
+ ),
+ # C4: the spheroid polar radius is negative.
+ # Expected: zero everywhere, warning that the polar radius must be
+ # positive.
+ (
+ "spheroidal_particle",
+ (numpy.array([1.0, 5.0]), 10.0, -5.0),
+ "'polar_radius' must be positive",
+ ),
+ # C5: the mean particle diameter of the lognormal distribution is
+ # zero.
+ # Expected: zero everywhere, warning that the diameter must be
+ # positive.
+ (
+ "lognormal_spherical_distribution",
+ (numpy.array([1.0, 2.0]), 0.0, 1.0),
+ "'particle_diameter' must be positive",
+ ),
+ # C6: the width of the lognormal distribution is negative. A negative
+ # width is meaningless rather than a limiting case, so it must not be
+ # quietly treated as a distribution of zero width.
+ # Expected: zero everywhere, warning that the width must not be
+ # negative.
+ (
+ "lognormal_spherical_distribution",
+ (numpy.array([1.0, 5.0, 10.0]), 10.0, -2.0),
+ "'particle_diameter_sigma' must not be negative",
+ ),
+ # C7: the mean particle diameter is too small for the requested
+ # distribution width, which puts the underlying lognormal mean below
+ # zero and invalidates the closed form.
+ # Expected: zero everywhere, warning that the diameter is too small
+ # for the width.
+ (
+ "lognormal_spherical_distribution",
+ (numpy.array([1.0, 5.0]), 2.0, 4.0),
+ "'particle_diameter' is too small for 'particle_diameter_sigma'",
+ ),
+ # C8: the sheet thickness is zero.
+ # Expected: zero everywhere, warning that the thickness must be
+ # positive.
+ (
+ "sheet_particle",
+ (numpy.array([1.0, 2.0]), 0.0),
+ "'sheet_thickness' must be positive",
+ ),
+ # C9: the sheet thickness is negative.
+ # Expected: zero everywhere, warning that the thickness must be
+ # positive.
+ (
+ "sheet_particle",
+ (numpy.array([1.0, 2.0]), -1.0),
+ "'sheet_thickness' must be positive",
+ ),
+ # C10: the shell thickness is zero, which would otherwise make the
+ # normalizing denominator vanish and return an unattenuated one at
+ # every r.
+ # Expected: zero everywhere, warning that the thickness must be
+ # positive.
+ (
+ "spherical_shell",
+ (numpy.array([1.0, 5.0, 10.0]), 10.0, 0.0),
+ "'thickness' must be positive",
+ ),
+ # C11: the shell thickness is negative, which would otherwise return
+ # a smoothly varying branch with values below negative one that an
+ # optimizer can mistake for a real solution.
+ # Expected: zero everywhere, warning that the thickness must be
+ # positive.
+ (
+ "spherical_shell",
+ (numpy.array([1.0, 5.0, 10.0]), 10.0, -5.0),
+ "'thickness' must be positive",
+ ),
+ # C12: the shell inner radius is negative.
+ # Expected: zero everywhere, warning that the radius must not be
+ # negative.
+ (
+ "spherical_shell",
+ (numpy.array([1.0, 5.0, 10.0]), -20.0, 5.0),
+ "'radius' must not be negative",
+ ),
+ ],
+)
+def test_non_physical_input_returns_zero_and_warns(
+ function_name,
+ input_args,
+ expected_requirement,
+ reset_characteristic_function_warnings,
+):
+ function = getattr(cf, function_name)
+ input_r = input_args[0]
+ expected_remedy = (
+ "The characteristic function was set to zero for every r, which "
+ "flattens the fit residual so a refinement cannot recover on its own. "
+ "Please use a physically meaningful starting value, or keep the "
+ "parameter in range with FitRecipe.add_soft_bounds."
+ )
+ expected_msg = (
+ f"In '{function_name}', {expected_requirement}. {expected_remedy}"
+ )
+ with pytest.warns(RuntimeWarning, match=re.escape(expected_msg)):
+ actual_characteristic_function = function(*input_args)
+ npt.assert_allclose(
+ actual_characteristic_function, numpy.zeros_like(input_r)
+ )
+
+
+@pytest.mark.parametrize(
+ "function_name, input_args, expected_characteristic_function",
+ [
+ # C1: a lognormal distribution of zero width is the sphere limit
+ # rather than a mistake.
+ # Expected: matches spherical_particle for the same diameter, with no
+ # warning.
+ (
+ "lognormal_spherical_distribution",
+ (numpy.array([1.0, 5.0, 10.0]), 10.0, 0.0),
+ cf.spherical_particle(numpy.array([1.0, 5.0, 10.0]), 10.0),
+ ),
+ # C2: a shell of zero inner radius is a solid sphere rather than a
+ # mistake.
+ # Expected: matches spherical_particle for twice the thickness, with
+ # no warning.
+ (
+ "spherical_shell",
+ (numpy.array([1.0, 5.0, 10.0]), 0.0, 5.0),
+ cf.spherical_particle(numpy.array([1.0, 5.0, 10.0]), 10.0),
+ ),
+ ],
+)
+def test_limiting_input_does_not_warn(
+ function_name,
+ input_args,
+ expected_characteristic_function,
+ reset_characteristic_function_warnings,
+):
+ function = getattr(cf, function_name)
+ with warnings.catch_warnings():
+ warnings.simplefilter("error")
+ actual_characteristic_function = function(*input_args)
+ npt.assert_allclose(
+ actual_characteristic_function, expected_characteristic_function
+ )
+
+
+def test_non_physical_input_warns_only_once(
+ reset_characteristic_function_warnings,
+):
+ # A refinement evaluates the characteristic function on every iteration,
+ # so a function that warns on each call would flood the user. The
+ # suppression must not rely on the warnings registry, which callers can
+ # disable, so this asserts under the "always" filter.
+ with warnings.catch_warnings(record=True) as raised_warnings:
+ warnings.simplefilter("always")
+ for _ in range(5):
+ cf.spherical_particle(numpy.array([1.0, 2.0]), -10.0)
+ actual_warning_count = len(raised_warnings)
+ expected_warning_count = 1
+ assert actual_warning_count == expected_warning_count
+
+
+if __name__ == "__main__":
+ unittest.main()