From 7217eb9601bc220d00a2af28c6fb3baf7592d6a8 Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Wed, 29 Jan 2025 15:11:39 -0500 Subject: [PATCH 1/4] Round to nearest integer for star polymer arms. Fixes #298 --- sasmodels/models/star_polymer.c | 19 +++++++++++-------- sasmodels/models/star_polymer.py | 4 ++-- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/sasmodels/models/star_polymer.c b/sasmodels/models/star_polymer.c index 13849789..2c687b6e 100644 --- a/sasmodels/models/star_polymer.c +++ b/sasmodels/models/star_polymer.c @@ -2,21 +2,21 @@ double form_volume(void); double Iq(double q, double rg_sq, double arms); -static double star_polymer_kernel(double q, double rg_sq, double arms) +static double star_polymer_kernel(double q, double rg_sq, int n_arms) { // Given the kernel function: // I = 2(v + (e^-v - 1) + (a-1)/2 (e^-v - 1)^2) / av^2 // At q = 0, u_2 and v are 0 so term1 = term2 = 0, hence 0/0 = nan - // We can clean up the values around low q using Taylor expansions for + // We can clean up the values at low q using Taylor expansions for // T1 = (e^-v - 1)/v = [-1 +1/2 -1/6 +1/24 -1/120 +1/720 ...] // T2 = [(e^-v - 1)/v]^2 = [+1 -1 +7/12 -1/4 +31/360 -1/40 ...] // Rewriting // I = 2/av + 2/av (e^-v - 1)/v + (a-1)/a ((e^v - 1)/v)^2 // Since the first term of T1 is -1 this cancels the leading 2/av. // Muliplying the remaining terms by 2/v yields T1' = 2*T1[1:] - // I = 1/a T1'(v) + (1 - 1/a) T2(v) + // I = 1/a T1'(v) + (1 - 1/a) T2(v) = (T1'(v) - T2(v))/a + T2(v) // T1' = [+1 -1/3 +1/12 -1/60 +1/360 -1/2520] - // Substituting v=0 we get lim[q->0] I = 1/a + (1 - 1/a) = 1 as expected. + // Substituting v=0 we get lim[q->0] I = (1 - 1)/a + 1 = 1 as expected. // After testing it turns out that term2 is numerically stable, so we // only use the Taylor series for T1', and protect T2 against v = 0. @@ -32,7 +32,7 @@ static double star_polymer_kernel(double q, double rg_sq, double arms) // Note: arms should be one or more const double u_sq = rg_sq * q * q; - const double v = u_sq * arms / (3.0 * arms - 2.0); + const double v = (u_sq * n_arms) / (3 * n_arms - 2); double term1; if (u_sq < _SPLIM) { term1 = 1. + v*(-1./3. + v*(1./12. + v*(-1./60. + v*(1./360. + v*(-1./2520))))); @@ -40,8 +40,7 @@ static double star_polymer_kernel(double q, double rg_sq, double arms) term1 = 2.0 * (v + expm1(-v)) / (v * v); } double term2 = square(expm1(-v)/v); // we trap q == 0 above, so v > 0 here - return term1/arms + term2*(1.0 - 1.0/arms); - + return (term1 - term2)/n_arms + term2; } double form_volume(void) @@ -51,5 +50,9 @@ double form_volume(void) double Iq(double q, double rg_sq, double arms) { - return star_polymer_kernel(q, rg_sq, arms); + // TODO: round or truncate? + // rounding: stacked_disk pearl_necklace multilayer_vesicle linear_pearls + // lamellar_stack_caille lamellar_hg_stack_caille + // truncation: lamellar_stack_paracrystal + return star_polymer_kernel(q, rg_sq, int(arms + 0.5)); } diff --git a/sasmodels/models/star_polymer.py b/sasmodels/models/star_polymer.py index e71be621..ab2f8f36 100644 --- a/sasmodels/models/star_polymer.py +++ b/sasmodels/models/star_polymer.py @@ -78,7 +78,7 @@ # pylint: disable=bad-whitespace, line-too-long # ["name", "units", default, [lower, upper], "type","description"], parameters = [["rg_squared", "Ang^2", 100.0, [0.0, inf], "", "Ensemble radius of gyration SQUARED of the full polymer"], - ["arms", "", 3, [1.0, 6.0], "", "Number of arms in the model"], + ["arms", "", 3, [1.0, 6.0], "", "Number of arms in the model (integer)"], ] # pylint: enable=bad-whitespace, line-too-long @@ -104,4 +104,4 @@ def random(): }, 1.0, 2.53575888234], ] # 2016-03-23 RKH edited docs, would this better use rg not rg^2 ? Numerical noise at extremely small q.rg -# 2025-01-17 PAK Fixed numerical noise. +# 2025-01-17 PAK Fixed numerical noise. Use nearest integer for number of arms. From 4945775aab9d0ae084fa1f5fe52b45f324fe7d68 Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Wed, 29 Jan 2025 15:26:50 -0500 Subject: [PATCH 2/4] pocl on CI pickier about C++ constructs in C code --- sasmodels/models/star_polymer.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sasmodels/models/star_polymer.c b/sasmodels/models/star_polymer.c index 2c687b6e..b819ad57 100644 --- a/sasmodels/models/star_polymer.c +++ b/sasmodels/models/star_polymer.c @@ -54,5 +54,5 @@ double Iq(double q, double rg_sq, double arms) // rounding: stacked_disk pearl_necklace multilayer_vesicle linear_pearls // lamellar_stack_caille lamellar_hg_stack_caille // truncation: lamellar_stack_paracrystal - return star_polymer_kernel(q, rg_sq, int(arms + 0.5)); + return star_polymer_kernel(q, rg_sq, (int)(arms + 0.5)); } From f8f77796a6640e067eb3f2aabd3514d849adc6ed Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Wed, 29 Jan 2025 15:32:38 -0500 Subject: [PATCH 3/4] test output for non-integer arms has now changed --- sasmodels/models/star_polymer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sasmodels/models/star_polymer.py b/sasmodels/models/star_polymer.py index ab2f8f36..be038204 100644 --- a/sasmodels/models/star_polymer.py +++ b/sasmodels/models/star_polymer.py @@ -96,7 +96,7 @@ def random(): tests = [[{'rg_squared': 2.0, 'arms': 3.3, - }, 0.5, 0.851646091108], + }, 0.5, 0.8518871750393033], [{'rg_squared': 1.0, 'arms': 2.0, From 509c642964858b9eda1326b3cbeb23a143b1173a Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Wed, 29 Jan 2025 15:40:26 -0500 Subject: [PATCH 4/4] test output for non-integer arms has now changed --- sasmodels/models/star_polymer.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/sasmodels/models/star_polymer.py b/sasmodels/models/star_polymer.py index be038204..8464e28b 100644 --- a/sasmodels/models/star_polymer.py +++ b/sasmodels/models/star_polymer.py @@ -95,7 +95,11 @@ def random(): return pars tests = [[{'rg_squared': 2.0, - 'arms': 3.3, + 'arms': 3.3, # should equal the value for arms=3 in the next test + }, 0.5, 0.8518871750393033], + + [{'rg_squared': 2.0, + 'arms': 3.0, }, 0.5, 0.8518871750393033], [{'rg_squared': 1.0,