Opened 6 years ago

Last modified 6 years ago

#1121 new enhancement

validation of beta(Q) and other calculations against external code

Reported by: richardh Owned by:
Priority: critical Milestone: sasmodels 1.0
Component: sasmodels Keywords:
Cc: Work Package: Beta Approximation Project

Description

Please document here detail results of testing against sasfit, Fish, own matlab codes etc, especially noting where issue were found and how these were resolved.

This will provide background information for sasview & sasmodels documentation.

Attachments (7)

sphere1.pdf (66.2 KB) - added by richardh 6 years ago.
sphere2.pdf (77.3 KB) - added by richardh 6 years ago.
sphere2a.pdf (60.3 KB) - added by richardh 6 years ago.
sphere3.pdf (73.8 KB) - added by richardh 6 years ago.
ellipsoid4..pdf (73.6 KB) - added by richardh 6 years ago.
ellipsoid5.pdf (48.0 KB) - added by richardh 6 years ago.
ellipsoid6bad.pdf (82.1 KB) - added by richardh 6 years ago.

Download all attachments as: .zip

Change History (11)

comment:1 Changed 6 years ago by pkienzle

sasmodels/explore/beta/sasfit_compare.py in the beta_approx branch contains code for computing <F> <F2> P S I=P*S S_eff and I_beta = P*S_eff in a way that is compatible with sasview 4.2, Yun's matlab code and sasfit.

As of this writing, the following alogithm is used, as cribbed from the ellipsoid_pe function:

# integrate over polydispersity in shape
F1 = F2 = 0
F_norm = 0
for wk, pk in distribution:
    # integrate over angles with u = cos(theta) substitution and z = 2*u - 1 
    # so that we have the equivalent of sum(F(t)*sin(t) dt) from 0 to pi/2 when
    # using gauss-legendre integration values and weights from -1 to 1.
    F1k = F2k = 0
    for w,z in gaussian_weights:
        form = contrast*volume*F(q, pk, cos_theta = (z+1)/2)
        F1k += w * form / 2
        F2k += w * form*form / 2

    # accumulate 1D patterns over polydispersity
    if sasfit:
        F1 += wk * F1k
        F2 += wk * F2k 
        F_norm += wk
    elif sasview or yun:
        F1 += wk * volume * F1k/volume
        F2 += wk * volume * F2k/volume
        F_norm += wk * volume

F1 = F1/F_norm
F2 = F2/F_norm
Sq = S(q)
beta = F1**2/F2

if sasfit:
    Sq_eff = 1 + beta * (Sq - 1)
    Pq = F2
elif sasview:
    Sq_eff = undefined
    Pq = F2*1e-4*volfraction
elif yun:
    Sq_eff = 1 + F_norm*beta * (Sq - 1)
    Pq = F2/volume

Iq = Pq*Sq
Iq_beta = Pq*Sq_eff

The match to sasfit is not very good (in the order of 1% relative error), but probably good enough that we know what calculation is being performed even though we are using different numeric procedures to perform it.

Still haven't seen Yun's results for gaussian dispersion of radii; may need to adapt the algorithm a bit once these are available.

Version 0, edited 6 years ago by pkienzle (next)

comment:2 Changed 6 years ago by richardh

For information, Greg's original fork of sasmodels is here:

https://github.com/gregsuczewski/sasmodels

comment:3 Changed 6 years ago by richardh

Some progress today attempting to validate the beta(Q) calcs in qt5 ESS_GUI branch of sasview and beta_approx branch of sasmodels

I have started with data files oddly called richard_test etc from sasmodels/explore/beta/data which can be loaded into sasview, code fragment below from sasfit_compare.py gives the details of their contents, polydisperse spheres and ellipsoid time hard sphere, with or without beta(Q) included.

def compare_sasfit_sphere_schulz():

#radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1
#We have scaled the output from sasfit by 1e-4*volume*volfraction
#0.10050378152592121
pars = {

'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz',
'sld': 4, 'sld_solvent': 1,
'volfraction': 0.3,

}

Q, IQ = load_sasfit(data_file('richard_test.txt'))
Q, IQSD = load_sasfit(data_file('richard_test2.txt'))
Q, IQBD = load_sasfit(data_file('richard_test3.txt'))


def compare_sasfit_ellipsoid_schulz():

#polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1
#Effective radius =13.1353356684
#We have scaled the output from sasfit by 1e-4*volume*volfraction
#0.10050378152592121
pars = {

'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz',
'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz',
'sld': 4, 'sld_solvent': 1,
'volfraction': 0.3, 'radius_effective': 13.1353356684,

}

Q, IQ = load_sasfit(data_file('richard_test4.txt'))
Q, IQSD = load_sasfit(data_file('richard_test5.txt'))
Q, IQBD = load_sasfit(data_file('richard_test6.txt'))

Thus far sasview appears to produce results consistent with these files from sasfit, bar issues with the overall scale which may require some thought.

I will attach some pdf reports to demonstrate, however note that the reports are alas not yet reporting polydispersity!

However in one case I got cranky results, see ellipsoid6bad where note the poor fit at small Q. Loading the data and setting up again in a new fit tab gave OK results see ellipsoid6better, this is rather worrying. Changing model in the original fit tab then back again to correct model seemed to sort the issue, which might be due to plotting problems?

On appending P(Q) to the plots in some cases, note that P(Q) is much smaller than you might expect due to "scale" being large, even if scale is set to 1.0 and sld adjusted the appended P(Q)is still nowhere near the fit. This needs some thought but may be a consequence of our single scale (or volume) parameter.

Changed 6 years ago by richardh

Changed 6 years ago by richardh

Changed 6 years ago by richardh

Changed 6 years ago by richardh

Changed 6 years ago by richardh

Changed 6 years ago by richardh

Changed 6 years ago by richardh

comment:4 Changed 6 years ago by pkienzle

  • Component changed from SasView to sasmodels
  • Milestone changed from SasView Next Release +1 to sasmodels 1.0
Note: See TracTickets for help on using tickets.