Skip to content

Single kernel call for volume and intensity #3963

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jellybean2004 wants to merge 3 commits into
base: SD_calc_reformat
Choose a base branch
from
Open

Single kernel call for volume and intensity #3963

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d921180
Single call for volume and intensity
jellybean2004 May 21, 2026
7adbbac
Reviewer's changes
jellybean2004 May 26, 2026
22beeea
Replace ellipse_volume with self._volumes
jellybean2004 May 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 30 additions & 28 deletions src/sas/sascalc/size_distribution/SizeDistribution.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@

from sasdata.dataloader.data_info import Data1D
from sasmodels.core import load_model
from sasmodels.direct_model import DirectModel
from sasmodels.direct_model import call_Fq

from sas.sascalc.size_distribution.maxEnt_method import maxEntMethod

Expand Down Expand Up @@ -114,6 +114,7 @@ def background_fit(
# Fit only scale
def fit_func(x: npt.ArrayLike, b: float) -> npt.ArrayLike:
return line_func(x, b, power)

init_guess = linearized_data.y[0]

else:
Expand All @@ -138,16 +139,6 @@ def fit_func(x: npt.ArrayLike, b: float) -> npt.ArrayLike:
return param_result, param_err


def ellipse_volume(rp: float, re: float) -> float:
"""
Calculate the volume of an ellipsoid given the polar and equatorial radii.
:param rp: polar radius
:param re: equatorial radius
:return: volume of the ellipsoid
"""
return (4.0 * np.pi / 3.0) * rp * re**2.0


class sizeDistribution:
"""Class for performing size distribution analysis using the Maximum Entropy method."""

Expand Down Expand Up @@ -177,6 +168,7 @@ def __init__(self, data: Data1D):
self._resolution: float | None = None

self.model_matrix: np.ndarray | None = None
self._volumes: np.ndarray | None = None

# Advanced parameters for MaxEnt
self._iterMax: int = 5000
Expand All @@ -188,7 +180,6 @@ def __init__(self, data: Data1D):

self._bin_edges: np.ndarray = np.array([])
self._binDiff: np.ndarray = np.array([])
self._volumes: np.ndarray | None = None

# Return Values after the MaxEnt fit
self.BinMagnitude_maxEnt: np.ndarray = np.array([], dtype=float)
Expand Down Expand Up @@ -483,37 +474,48 @@ def generate_model_matrix(self, moddata: Data1D) -> None:
:param moddata: Data1D object that has the data trimmed depending on background
subtraction or power law subtracted from the data. Also self.qMin and self.qMax.
"""
model = load_model(self.model)

pars = {
"sld": self.contrast,
"sld_solvent": 0.0,
"background": 0.0,
"scale": 1.0,
}

kernel = DirectModel(moddata, model)

intensities = []
volumes = []

# Build a single Kernel to compute both intensity and volume per bin
model_obj = load_model(self.model)
calculator = model_obj.make_kernel((moddata.x,))

for bin in self.bins:
pars["radius_equatorial"] = bin
pars["radius_polar"] = bin * self.aspectRatio
intensities.append(kernel(**pars))
p = pars.copy()
p["radius_equatorial"] = bin
p["radius_polar"] = bin * self.aspectRatio

_, Fsq, _, volume, volume_ratio = call_Fq(calculator, p)

# Compute intensity using kernel convention: combined_scale = scale / shell_volume
Copy link
Copy Markdown
Contributor

@pkienzle pkienzle May 22, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is duplicating code from within sasmodels.kernel.Kernel.Iq. This isn't avoidable in the current implementation, but we should consider this use case when addressing sasmodels #182.

The idea is to return intermediate and derived values from the calculation back to the caller. The current hack to get $P(Q)$ and $S(Q)$ from P@S is that sasmodels.product sets up a kernel.results() function that returns a dictionary with keys P(Q), S(Q), volume, volume_ratio and radius_effective after calling kernel.Iq(). We could extend this to the base Kernel class in sasmodels. Then using call_kernel instead of call_Fq we can get the volume and volume ratio from calculator.results().

scale_val = p.get("scale", 1.0)
background_val = p.get("background", 0.0)
combined_scale = scale_val / volume

intensity = combined_scale * Fsq + background_val
intensities.append(intensity)

# Volume ratio is 1 for solid particles
volume *= volume_ratio
volumes.append(volume)

self.model_matrix = np.vstack(intensities).T
self._volumes = np.array(volumes)

def calc_volume_weighted_dist(self, binmag: np.ndarray) -> None:
"""
This is not used right now.
Calculate the volume weighted distribution.
"""
if self.logbin:
radbins = np.logspace(np.log10(self.diamMin), np.log10(self.diamMax), self.nbins + 1, True) * 0.5

else:
radbins = np.linspace(self.diamMin, self.diamMax, self.nbins + 1, True) * 0.5

self.volume_bins = ellipse_volume(self.aspectRatio * radbins, radbins)
self.volume_bins = self._volumes
self.vbin_diff = np.diff(self.volume_bins)
self.volume_bins = self.volume_bins[:-1] + self.vbin_diff * 0.5
self.volume_fraction = binmag * self.volume_bins / (2.0 * self.vbin_diff)
Expand Down Expand Up @@ -669,12 +671,12 @@ def calculate_statistics(self, bin_mag: npt.ArrayLike) -> None:
"""
Calculate statistics from the MaxEnt results, including volume fraction cumulative distribution function (CDF),
number distribution, and related statistics such as mean, median, and mode.

:param bin_mag: list of bin magnitudes from the MaxEnt fits
"""
bin_mag = np.asarray(bin_mag)
maxent_cdf_array = integrate.cumulative_trapezoid(bin_mag / (2.0 * self._binDiff), 2.0 * self.bins, axis=1)
self.BinMag_numberDist = self.BinMagnitude_maxEnt / ellipse_volume(self.aspectRatio * self.bins, self.bins)
self.BinMag_numberDist = self.BinMagnitude_maxEnt / self._volumes

rvdist = stats.rv_histogram(
(self.BinMagnitude_maxEnt, self._bin_edges * 2.0), density=True
Expand Down
Loading