From 8a160944386a795605b7852450bc48c07c353d4e Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Fri, 1 May 2026 12:12:48 -0400 Subject: [PATCH] verify accuracy of fractal S(q) for small q --- explore/precision.py | 43 +++++++++++++++++++++++++++++++ sasmodels/models/fractal.py | 2 +- sasmodels/models/lib/fractal_sq.c | 24 +++++++---------- sasmodels/special/__init__.py | 21 ++++++++++++++- 4 files changed, 73 insertions(+), 17 deletions(-) diff --git a/explore/precision.py b/explore/precision.py index fdc009b78..6951e38c2 100755 --- a/explore/precision.py +++ b/explore/precision.py @@ -183,6 +183,7 @@ def compare(self, x, precision, target, linear=False, diff="relative"): pylab.semilogx(x, target, '-', label="true value") if linear: pylab.xscale('linear') + #pylab.yscale('linear') def plotdiff(x, target, actual, label, diff): """ @@ -192,10 +193,20 @@ def plotdiff(x, target, actual, label, diff): """ if diff == "relative": err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd') + if (err == 0).any(): + if (err > 0).any(): + err[err == 0] = err[err > 0].min()/2 + else: + err[err == 0] = 1e-20 #err = np.clip(err, 0, 1) pylab.loglog(x, err, '-', label=label, alpha=0.7) elif diff == "absolute": err = np.array([abs(t-a) for t, a in zip(target, actual)], 'd') + if (err == 0).any(): + if (err > 0).any(): + err[err == 0] = err[err > 0].min()/2 + else: + err[err == 0] = 1e-20 pylab.loglog(x, err, '-', label=label, alpha=0.7) else: limits = np.min(target), np.max(target) @@ -608,6 +619,37 @@ def np_star_polymer(x, arms=3): xaxis="$(Q R_g)^2$ (unitless)", ) +from sasmodels.special import fractal_sq + + +def mp_fractal_sq(q, radius, fractal_dim, cor_length): + # Eq 16 from Texiera (1988), DOI:10.1107/S0021889888000263 + D = fractal_dim + if D == 0: + term = 1 + elif D == 1: + term = 1 if q == 0 else mp.atan(q*cor_length)/(q*radius) + elif q == 0: + term = mp.power(cor_length/radius, D)*mp.gamma(D+1) + else: + t1 = mp.power(q*radius, -D) + t2 = D*mp.gamma(D-1) / mp.power(1 + 1/(q*cor_length)**2, (D-1)/2) + t3 = mp.sin((D-1)*mp.atan(q*cor_length)) + term = t1*t2*t3 + + return 1 + term + +FRACTAL_RADIUS=30 +FRACTAL_DIM=1.6 +FRACTAL_CORR=421 +add_function( + name="fractal", + mp_function=lambda x: mp_fractal_sq(x, FRACTAL_RADIUS, FRACTAL_DIM, FRACTAL_CORR), + np_function=lambda x: fractal_sq(x, FRACTAL_RADIUS, FRACTAL_DIM, FRACTAL_CORR), + ocl_function=make_ocl(f"return fractal_sq(q, {FRACTAL_RADIUS}, {FRACTAL_DIM}, {FRACTAL_CORR});", + "fractal", source=["lib/sas_gamma.c", "lib/fractal_sq.c"]), +) + RADIUS=3000 LENGTH=30 @@ -686,6 +728,7 @@ def np_cyl(x): ocl_function=make_ocl(lanczos_gamma, "lgamma"), ) + replacement_expm1 = """\ double x = (double)q; // go back to float for single precision kernels // Adapted from the cephes math library. diff --git a/sasmodels/models/fractal.py b/sasmodels/models/fractal.py index c3dd60d8d..3544ff50d 100644 --- a/sasmodels/models/fractal.py +++ b/sasmodels/models/fractal.py @@ -43,7 +43,7 @@ References ---------- -#. J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 +#. J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 DOI:10.1107/S0021889888000263 Authorship and Verification ---------------------------- diff --git a/sasmodels/models/lib/fractal_sq.c b/sasmodels/models/lib/fractal_sq.c index 689c83bbe..c8b0de198 100644 --- a/sasmodels/models/lib/fractal_sq.c +++ b/sasmodels/models/lib/fractal_sq.c @@ -1,28 +1,22 @@ static double fractal_sq(double q, double radius, double fractal_dim, double cor_length) { - //calculate S(q), using Teixeira, Eq(15) - // mathematica query to check limiting conditions: - // lim x->0 of [ x gamma(x-1) sin(arctan(q c (x-1))) (q r)^(-x) (1 + 1/(q c)^2)^((1-x)/2) ] - // Note: gamma(x) may be unreliable for x<0, so the gamma(D-1) is risky. - // We instead transform D*gamma(D-1) into gamma(D+1)/(D-1). + // Calculate S(q) using Eq(16) from Teixeira (1988), DOI:10.1107/S0021889888000263 + // with terms rearranged so that limiting cases are easier to calculate. + // Values checked against 500 bit floating point using the original equations. double term; if (q == 0.) { const double D = fractal_dim; - term = pow(cor_length/radius, D)*sas_gamma(D+1.); + term = pow(cor_length/radius, fractal_dim)*sas_gamma(fractal_dim+1.); } else if (fractal_dim == 0.) { - term = 1.0; + term = 1.; } else if (fractal_dim == 1.) { term = atan(q*cor_length)/(q*radius); } else { - // q>0, D>0 - const double D = fractal_dim; - const double Dm1 = fractal_dim - 1.0; - // Note: for large Dm1, sin(Dm1*atan(q*cor_length) can go negative - const double t1 = sas_gamma(D+1.)/Dm1*sin(Dm1*atan(q*cor_length)); - const double t2 = pow(q*radius, -D); - const double t3 = pow(1.0 + 1.0/square(q*cor_length), -0.5*Dm1); + const double t1 = pow(cor_length/radius, fractal_dim) * sas_gamma(fractal_dim+1.); + const double t2 = sin((fractal_dim-1.)*atan(q*cor_length))/((fractal_dim-1)*q*cor_length); + const double t3 = pow(1. + square(q*cor_length), -0.5*(fractal_dim-1.)); term = t1 * t2 * t3; } - return 1.0 + term; + return 1. + term; } diff --git a/sasmodels/special/__init__.py b/sasmodels/special/__init__.py index 390173eed..e6a0ec8f5 100644 --- a/sasmodels/special/__init__.py +++ b/sasmodels/special/__init__.py @@ -14,7 +14,8 @@ # Functions to add to our standard set # C99 standard math library functions -from numpy import cos, inf, nan, sin +from numpy import acos, asin, atan, cos, inf, nan, pow, sin, tan +from scipy.special import gamma as sas_gamma from scipy.special import j1 as sas_J1 # erf, erfc, tgamma, lgamma **do not use** @@ -98,6 +99,24 @@ def sas_2J1x_x(x): retvalue[x == 0] = 1. return retvalue +def fractal_sq(q, radius, fractal_dim, cor_length): + # Calculate S(q) using Eq(16) from Teixeira (1988), DOI:10.1107/S0021889888000263 + # with terms rearranged so that limiting cases are easier to calculate. + # Values checked against 500 bit floating point using the original equations. + with np.errstate(all='ignore'): + if fractal_dim == 0: + term = np.ones_like(q) + elif fractal_dim == 1: + term = atan(q*cor_length)/(q*radius) + term[q == 0] = 1 + else: + t1 = pow(cor_length/radius, fractal_dim) * sas_gamma(fractal_dim+1) + t2 = sin((fractal_dim-1)*atan(q*cor_length))/((fractal_dim-1)*q*cor_length) + t2[q == 0] = 1 + t3 = pow(1 + square(q*cor_length), -(fractal_dim-1)/2) + term = t1*t2*t3 + + return 1 + term # Gaussians class Gauss: