From e4b2ff4349b18544c545ae83da79427334ccee35 Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Wed, 21 Oct 2026 16:39:53 -0400 Subject: [PATCH 1/6] rename functions, add DANSE license --- LICENSE_DANSE.rst | 34 ++ news/fix-cfs.rst | 24 ++ src/diffpy/cmipdf/characteristicfunctions.py | 420 ++++++++++++------- 3 files changed, 323 insertions(+), 155 deletions(-) create mode 100644 LICENSE_DANSE.rst create mode 100644 news/fix-cfs.rst diff --git a/LICENSE_DANSE.rst b/LICENSE_DANSE.rst new file mode 100644 index 0000000..c92ca2d --- /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 sb2896@columbia.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/news/fix-cfs.rst b/news/fix-cfs.rst new file mode 100644 index 0000000..e4c0382 --- /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:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/src/diffpy/cmipdf/characteristicfunctions.py b/src/diffpy/cmipdf/characteristicfunctions.py index bebf880..1755b4c 100644 --- a/src/diffpy/cmipdf/characteristicfunctions.py +++ b/src/diffpy/cmipdf/characteristicfunctions.py @@ -34,6 +34,8 @@ "SASCF", ] +import warnings + import numpy from numpy import arctan as atan from numpy import arctanh as atanh @@ -44,91 +46,133 @@ 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 +186,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 +201,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 +222,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 +242,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 +250,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. - - 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. +def lognormal_spherical_distribution( + r, particle_diameter, particle_diameter_sigma +): + """Compute the spherical nanoparticle characteristic function with + lognormal size distribution. - - 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 +268,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((particle_diameter_sigma/particle_diameter)^2 + 1) + mu = log(particle_diameter) - s^2/2 - s^2 = log((p_sigma/p_diameter)^2 + 1) - mu = log(p_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 @@ -390,10 +489,10 @@ def __init__(self, name, model): Parameters ---------- - name - A name for the SASCF - model - SASModel object this adapts. + name : str + The name for the SASCF. + model : BaseModel + The sas.models.BaseModel object this adapts. """ Calculator.__init__(self, name) @@ -418,7 +517,6 @@ def __init__(self, name, model): 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 @@ -428,8 +526,8 @@ def __call__(self, r): # # 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%. + # 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. @@ -485,6 +583,18 @@ def __call__(self, r): def erfc(x): + """Compute the complementary error function. + + Parameters + ---------- + x : array_like + The input value. + + Returns + ------- + numpy.ndarray + The complementary error function evaluated at `x`. + """ return 1.0 - erf(x) From bdf67bd4303719a5f06759c8ca1e4e2ae7f65b98 Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Wed, 21 Oct 2026 16:43:25 -0400 Subject: [PATCH 2/6] copy over cf testing --- tests/conftest.py | 35 ++ tests/test_characteristicfunctions.py | 581 ++++++++++++++++++++++++++ 2 files changed, 616 insertions(+) create mode 100644 tests/test_characteristicfunctions.py diff --git a/tests/conftest.py b/tests/conftest.py index e4544a1..7180881 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,5 +1,6 @@ import importlib.resources import json +from functools import lru_cache from pathlib import Path import pytest @@ -7,6 +8,40 @@ from diffpy.structure import PDFFitStructure +@lru_cache() +def has_sas(): + try: + import sas + import sasmodels + + del sas + del sasmodels + return True + except ImportError: + return False + + +@pytest.fixture(scope="session") +def sas_available(): + return has_sas() + + +@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..15985ea --- /dev/null +++ b/tests/test_characteristicfunctions.py @@ -0,0 +1,581 @@ +"""Tests for the nanoparticle characteristic functions and the SAS +calculator wrapper.""" + +# 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 + +try: + from sasmodels.sasview_model import find_model, load_standard_models +except ModuleNotFoundError: + pass + +import diffpy.cmipdf.characteristicfunctions as cf + +# # Global variables to be assigned in setUp +# cf = None +# Fixme: remove this code if imports don't break on testing + +# ---------------------------------------------------------------------------- + + +def testSphere(sas_available): + # if not sas_available: + # pytest.skip("sas package not available") + pytest.skip( + "sas characteristic functions not currently working, " + + "remove skip when our code is refactored to use the " + + "latest sasview API" + ) + + load_standard_models() + radius = 25 + # Calculate sphere cf from SphereModel + SphereModel = find_model("sphere") + model = SphereModel() + model.setParam("radius", radius) + ff = cf.SASCF("sphere", model) + r = numpy.arange(1, 60, 0.1, dtype=float) + fr1 = ff(r) + + # Calculate sphere cf analytically + fr2 = cf.spherical_particle(r, 2 * radius) + diff = fr1 - fr2 + res = numpy.dot(diff, diff) + res /= numpy.dot(fr2, fr2) + assert res == pytest.approx(0, abs=1e-4) + return + + +def testSpheroid(sas_available): + # if not sas_available: + # pytest.skip("sas package not available") + pytest.skip( + "sas characteristic functions not currently working, " + + "remove skip when our code is refactored to use the " + + "latest sasview API" + ) + + load_standard_models() + prad = 20.9 + erad = 33.114 + # Calculate cf from EllipsoidModel + EllipsoidModel = find_model("ellipsoid") + model = EllipsoidModel() + model.setParam("radius_polar", prad) + model.setParam("radius_equatorial", erad) + ff = cf.SASCF("spheroid", model) + r = numpy.arange(0, 100, 1 / numpy.pi, dtype=float) + fr1 = ff(r) + + # Calculate cf analytically + fr2 = cf.spheroidal_particle(r, erad, prad) + diff = fr1 - fr2 + res = numpy.dot(diff, diff) + res /= numpy.dot(fr2, fr2) + assert res == pytest.approx(0, abs=1e-4) + return + + +def testShell(sas_available): + # if not sas_available: + # pytest.skip("sas package not available") + pytest.skip( + "sas characteristic functions not currently working, " + + "remove skip when our code is refactored to use the " + + "latest sasview API" + ) + + load_standard_models() + radius = 19.2 + thickness = 7.8 + # Calculate cf from VesicleModel + VesicleModel = find_model("vesicle") + model = VesicleModel() + model.setParam("radius", radius) + model.setParam("thickness", thickness) + ff = cf.SASCF("vesicle", model) + r = numpy.arange(0, 99.45, 0.1, dtype=float) + fr1 = ff(r) + + # Calculate sphere cf analytically + fr2 = cf.spherical_shell(r, radius, thickness) + diff = fr1 - fr2 + res = numpy.dot(diff, diff) + res /= numpy.dot(fr2, fr2) + assert res == pytest.approx(0, abs=1e-4) + return + + +def testCylinder(sas_available): + # if not sas_available: + # pytest.skip("sas package not available") + pytest.skip( + "sas characteristic functions not currently working, " + + "remove skip when our code is refactored to use the " + + "latest sasview API" + ) + load_standard_models() + # Make sure cylinder works over different r-ranges. + radius = 100 + length = 30 + + CylinderModel = find_model("cylinder") + model = CylinderModel() + model.setParam("radius", radius) + model.setParam("length", length) + + ff = cf.SASCF("cylinder", model) + + r1 = numpy.arange(0, 10, 0.1, dtype=float) + r2 = numpy.arange(0, 50, 0.1, dtype=float) + r3 = numpy.arange(0, 100, 0.1, dtype=float) + r4 = numpy.arange(0, 500, 0.1, dtype=float) + + fr1 = ff(r1) + fr2 = ff(r2) + fr3 = ff(r3) + fr4 = ff(r4) + + d = fr1 - numpy.interp(r1, r2, fr2) + res12 = numpy.dot(d, d) + res12 /= numpy.dot(fr1, fr1) + assert res12 == pytest.approx(0, abs=1e-4) + + d = fr1 - numpy.interp(r1, r3, fr3) + res13 = numpy.dot(d, d) + res13 /= numpy.dot(fr1, fr1) + assert res13 == pytest.approx(0, abs=1e-4) + + d = fr1 - numpy.interp(r1, r4, fr4) + res14 = numpy.dot(d, d) + res14 /= numpy.dot(fr1, fr1) + assert res14 == pytest.approx(0, abs=1e-4) + + d = fr2 - numpy.interp(r2, r3, fr3) + res23 = numpy.dot(d, d) + res23 /= numpy.dot(fr2, fr2) + assert res23 == pytest.approx(0, abs=1e-4) + + d = fr2 - numpy.interp(r2, r4, fr4) + res24 = numpy.dot(d, d) + res24 /= numpy.dot(fr2, fr2) + assert res24 == pytest.approx(0, abs=1e-4) + + d = fr3 - numpy.interp(r3, r4, fr4) + res34 = numpy.dot(d, d) + res34 /= numpy.dot(fr3, fr3) + assert res34 == pytest.approx(0, abs=1e-4) + return + + +# End of class TestSASCF + +# ---------------------------------------------------------------------------- +# 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() From 0658ebcc81c5fb74e26845584681caf07effc765 Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Wed, 21 Oct 2026 16:44:48 -0400 Subject: [PATCH 3/6] pre-commit --- CHANGELOG.rst | 1 - src/diffpy/cmipdf/characteristicfunctions.py | 1 - 2 files changed, 2 deletions(-) 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/src/diffpy/cmipdf/characteristicfunctions.py b/src/diffpy/cmipdf/characteristicfunctions.py index 1755b4c..90575a0 100644 --- a/src/diffpy/cmipdf/characteristicfunctions.py +++ b/src/diffpy/cmipdf/characteristicfunctions.py @@ -45,7 +45,6 @@ from diffpy.srfit.fitbase.calculator import Calculator - # 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() From ee0a74bf93a242e9cbf5e63d61bdafe259480410 Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Wed, 21 Oct 2026 16:49:39 -0400 Subject: [PATCH 4/6] change simons email from columbia to ucsb --- CODE-OF-CONDUCT.rst | 2 +- LICENSE_DANSE.rst | 2 +- README.rst | 2 +- cookiecutter.json | 2 +- pyproject.toml | 4 ++-- 5 files changed, 6 insertions(+), 6 deletions(-) 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 index c92ca2d..9e65c74 100644 --- a/LICENSE_DANSE.rst +++ b/LICENSE_DANSE.rst @@ -6,7 +6,7 @@ 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 sb2896@columbia.edu. +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: 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/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'] From c5248bae11b539b041e28c24cd0be89b9672c99c Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Fri, 23 Oct 2026 16:15:03 -0400 Subject: [PATCH 5/6] test: remove SAS tests --- tests/conftest.py | 19 --- tests/test_characteristicfunctions.py | 167 +------------------------- 2 files changed, 1 insertion(+), 185 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index 7180881..0227883 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,6 +1,5 @@ import importlib.resources import json -from functools import lru_cache from pathlib import Path import pytest @@ -8,24 +7,6 @@ from diffpy.structure import PDFFitStructure -@lru_cache() -def has_sas(): - try: - import sas - import sasmodels - - del sas - del sasmodels - return True - except ImportError: - return False - - -@pytest.fixture(scope="session") -def sas_available(): - return has_sas() - - @pytest.fixture def reset_characteristic_function_warnings(): """Clear the record of already-issued non-physical input warnings. diff --git a/tests/test_characteristicfunctions.py b/tests/test_characteristicfunctions.py index 15985ea..c194194 100644 --- a/tests/test_characteristicfunctions.py +++ b/tests/test_characteristicfunctions.py @@ -1,5 +1,4 @@ -"""Tests for the nanoparticle characteristic functions and the SAS -calculator wrapper.""" +"""Tests for the nanoparticle characteristic functions.""" # FIXME: remove this line when `docformatter` fixes the blank line bug @@ -12,172 +11,8 @@ import numpy.testing as npt import pytest -try: - from sasmodels.sasview_model import find_model, load_standard_models -except ModuleNotFoundError: - pass - import diffpy.cmipdf.characteristicfunctions as cf -# # Global variables to be assigned in setUp -# cf = None -# Fixme: remove this code if imports don't break on testing - -# ---------------------------------------------------------------------------- - - -def testSphere(sas_available): - # if not sas_available: - # pytest.skip("sas package not available") - pytest.skip( - "sas characteristic functions not currently working, " - + "remove skip when our code is refactored to use the " - + "latest sasview API" - ) - - load_standard_models() - radius = 25 - # Calculate sphere cf from SphereModel - SphereModel = find_model("sphere") - model = SphereModel() - model.setParam("radius", radius) - ff = cf.SASCF("sphere", model) - r = numpy.arange(1, 60, 0.1, dtype=float) - fr1 = ff(r) - - # Calculate sphere cf analytically - fr2 = cf.spherical_particle(r, 2 * radius) - diff = fr1 - fr2 - res = numpy.dot(diff, diff) - res /= numpy.dot(fr2, fr2) - assert res == pytest.approx(0, abs=1e-4) - return - - -def testSpheroid(sas_available): - # if not sas_available: - # pytest.skip("sas package not available") - pytest.skip( - "sas characteristic functions not currently working, " - + "remove skip when our code is refactored to use the " - + "latest sasview API" - ) - - load_standard_models() - prad = 20.9 - erad = 33.114 - # Calculate cf from EllipsoidModel - EllipsoidModel = find_model("ellipsoid") - model = EllipsoidModel() - model.setParam("radius_polar", prad) - model.setParam("radius_equatorial", erad) - ff = cf.SASCF("spheroid", model) - r = numpy.arange(0, 100, 1 / numpy.pi, dtype=float) - fr1 = ff(r) - - # Calculate cf analytically - fr2 = cf.spheroidal_particle(r, erad, prad) - diff = fr1 - fr2 - res = numpy.dot(diff, diff) - res /= numpy.dot(fr2, fr2) - assert res == pytest.approx(0, abs=1e-4) - return - - -def testShell(sas_available): - # if not sas_available: - # pytest.skip("sas package not available") - pytest.skip( - "sas characteristic functions not currently working, " - + "remove skip when our code is refactored to use the " - + "latest sasview API" - ) - - load_standard_models() - radius = 19.2 - thickness = 7.8 - # Calculate cf from VesicleModel - VesicleModel = find_model("vesicle") - model = VesicleModel() - model.setParam("radius", radius) - model.setParam("thickness", thickness) - ff = cf.SASCF("vesicle", model) - r = numpy.arange(0, 99.45, 0.1, dtype=float) - fr1 = ff(r) - - # Calculate sphere cf analytically - fr2 = cf.spherical_shell(r, radius, thickness) - diff = fr1 - fr2 - res = numpy.dot(diff, diff) - res /= numpy.dot(fr2, fr2) - assert res == pytest.approx(0, abs=1e-4) - return - - -def testCylinder(sas_available): - # if not sas_available: - # pytest.skip("sas package not available") - pytest.skip( - "sas characteristic functions not currently working, " - + "remove skip when our code is refactored to use the " - + "latest sasview API" - ) - load_standard_models() - # Make sure cylinder works over different r-ranges. - radius = 100 - length = 30 - - CylinderModel = find_model("cylinder") - model = CylinderModel() - model.setParam("radius", radius) - model.setParam("length", length) - - ff = cf.SASCF("cylinder", model) - - r1 = numpy.arange(0, 10, 0.1, dtype=float) - r2 = numpy.arange(0, 50, 0.1, dtype=float) - r3 = numpy.arange(0, 100, 0.1, dtype=float) - r4 = numpy.arange(0, 500, 0.1, dtype=float) - - fr1 = ff(r1) - fr2 = ff(r2) - fr3 = ff(r3) - fr4 = ff(r4) - - d = fr1 - numpy.interp(r1, r2, fr2) - res12 = numpy.dot(d, d) - res12 /= numpy.dot(fr1, fr1) - assert res12 == pytest.approx(0, abs=1e-4) - - d = fr1 - numpy.interp(r1, r3, fr3) - res13 = numpy.dot(d, d) - res13 /= numpy.dot(fr1, fr1) - assert res13 == pytest.approx(0, abs=1e-4) - - d = fr1 - numpy.interp(r1, r4, fr4) - res14 = numpy.dot(d, d) - res14 /= numpy.dot(fr1, fr1) - assert res14 == pytest.approx(0, abs=1e-4) - - d = fr2 - numpy.interp(r2, r3, fr3) - res23 = numpy.dot(d, d) - res23 /= numpy.dot(fr2, fr2) - assert res23 == pytest.approx(0, abs=1e-4) - - d = fr2 - numpy.interp(r2, r4, fr4) - res24 = numpy.dot(d, d) - res24 /= numpy.dot(fr2, fr2) - assert res24 == pytest.approx(0, abs=1e-4) - - d = fr3 - numpy.interp(r3, r4, fr4) - res34 = numpy.dot(d, d) - res34 /= numpy.dot(fr3, fr3) - assert res34 == pytest.approx(0, abs=1e-4) - return - - -# End of class TestSASCF - # ---------------------------------------------------------------------------- # 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 From 93387c53b2319d237ab9d4b83646e625e45cdf0e Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Fri, 23 Oct 2026 16:19:09 -0400 Subject: [PATCH 6/6] remove SAS characteristic functions --- news/fix-cfs.rst | 2 +- src/diffpy/cmipdf/characteristicfunctions.py | 124 +------------------ 2 files changed, 2 insertions(+), 124 deletions(-) diff --git a/news/fix-cfs.rst b/news/fix-cfs.rst index e4c0382..4c23e04 100644 --- a/news/fix-cfs.rst +++ b/news/fix-cfs.rst @@ -13,7 +13,7 @@ **Removed:** -* +* Remove the non-functional ``SASCF`` calculator from ``characteristicfunctions``. SAS characteristic functions now live in ``diffpy.cmisas``. **Fixed:** diff --git a/src/diffpy/cmipdf/characteristicfunctions.py b/src/diffpy/cmipdf/characteristicfunctions.py index 90575a0..90eab31 100644 --- a/src/diffpy/cmipdf/characteristicfunctions.py +++ b/src/diffpy/cmipdf/characteristicfunctions.py @@ -31,7 +31,6 @@ "lognormal_spherical_distribution", "sheet_particle", "spherical_shell", - "SASCF", ] import warnings @@ -39,12 +38,9 @@ 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 - # 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() @@ -463,124 +459,6 @@ def spherical_shell(r, radius, thickness): 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. - - Attributes - ---------- - _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. - - Parameters - ---------- - name : str - The name for the SASCF. - model : BaseModel - The sas.models.BaseModel 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): """Compute the complementary error function.