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
total_weight = total_volume = 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
    F1 += wk * F1k
    F2 += wk * F2k 
    total_weight += wk
    total_volume += wk * volume

F1 = F1/total_weight
F2 = F2/total_weight
average_volume = total_volume/total_weight

if sasfit:
    Pq = F2
elif sasview:
    Pq = F2/average_volume*1e-4*volfraction
elif yun:
    Pq = F2/average_volume*1e8*volfraction  # slds not scaled by 1e-6

Sq = S(q)
beta = F1**2/F2
Sq_eff = 1 + beta * (Sq - 1)
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.

Last edited 6 years ago by pkienzle (previous) (diff)

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.