SasView / sasmodels

Package for calculation of small angle scattering models using OpenCL.
BSD 3-Clause "New" or "Revised" License
15 stars 27 forks source link

Show scaled distribution on weights plot #556

Closed pkienzle closed 1 year ago

pkienzle commented 1 year ago

The sasmodels -weights plot is not showing up and is not weighted correctly.

Currently it is plotting weight associated with each polydispersity point in the distribution $P$ as $w_k = P(x_k) / \sum P(x_k)$ rather than plotting the distribution itself.

Assuming $P(x)$ is approximately constant over the interval, $w_k$ is the weight associated with the interval from $k$ to $k+1$ with $wk \approx \int{x{k}'}^{x{k+1}'} P(x)\ {\rm d}x$ in the limit of large $n$. That is, $w_k$ is the area of a rectangle of height $P(x) \approx w_k / Δx_k$.

We are using the central point of equally spaced intervals $Δx$ so that a symmetric distribution has its peak at the middle $x$ value when there is an odd number of points. This means we have $n+1$ endpoints with $x_k' = xk - ½Δx$ for $k<=n$ and $x{n+1}' = x_n + ½Δx$ for the distribution $({\bf x}, {\bf w})$.

So weights.py should set dx = x[1]-x[0] and use plot(x,w/dx) when plotting the distribution.

This will automatically handle changes in the distribution width and the number of points, as well as cutoff below zero. It does not handle distribution cutoff with the model definition, such as the barbell model with radius_bell >= radius.

pkienzle commented 1 year ago

After the above change to the plotted value it shows the following.

Half the width leads to twice the height:

$ sascomp ellipsoid -poly -noplot -weights -profile radius_equatorial=radius_polar radius_polar_pd=0.1 radius_equatorial_pd=0.2

image

Changing n gives the same probability:

$ sascomp ellipsoid -poly -pars -noplot -weights -profile radius_equatorial=radius_polar radius_polar_pd=0.1 radius_equatorial_pd=0.1 radius_polar_pd_n=15

image

It even adjusts the height for cutoff below zero:

$ sascomp ellipsoid -poly -pars -noplot -weights -profile radius_equatorial=radius_polar/5 radius_polar_pd=0.2 radius_equatorial_pd=1.0

image
pkienzle commented 1 year ago

Since the option is -weights and that is what we are currently plotting we may want to leave things as they are. The weights are more useful than the probabilities within the context of sascomp, where the primary purpose is debugging numerical issues and developing new models.

I will therefore close this issue. SasView itself has a different audience and may want to plot distribution rather than weights.

It would be useful to label the y axis with the value being plotted. The x axis is harder to label since it may be mixi angular and length dispersity and the lengths may have different units. It is not worth addressing this complexity for a development tool.