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

[RichardHeenan]

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.

Migrated from http://trac.sasview.org/ticket/1121

{
    "status": "new",
    "changetime": "2018-10-30T22:43:13",
    "_ts": "2018-10-30 22:43:13.383456+00:00",
    "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.\n\nThis will provide background information for sasview & sasmodels documentation.",
    "reporter": "richardh",
    "cc": "",
    "resolution": "",
    "workpackage": "Beta Approximation Project",
    "time": "2018-07-01T21:39:33",
    "component": "sasmodels",
    "summary": "validation of beta(Q) and other calculations against external code",
    "priority": "critical",
    "keywords": "",
    "milestone": "sasmodels 1.0",
    "owner": "",
    "type": "enhancement"
}

[pkienzle]

Trac update at 2018/07/05 16:57:16:

• pkienzle changed _comment0 from:

sasmodels/explore/beta/sasfit_compare.py in the beta_approx branch contains code for computing <F^2^> P S I=PS S_eff and I_beta = PS_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 = 2u - 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 = contrastvolumeF(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 = F21e-4volfraction
elif yun:
Sq_eff = 1 + F_norm*beta * (Sq - 1)
Pq = F2/volume

Iq = PqSq
Iq_beta = PqSq_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.

to:

1530833592125933

• pkienzle changed _comment1 from:

sasmodels/explore/beta/sasfit_compare.py in the beta_approx branch contains code for computing <F^2^> P S I=PS S_eff and I_beta = PS_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 = 2u - 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 = contrastvolumeF(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
Sq = S(q)
beta = F1**2/F2
Sq_eff = 1 + beta * (Sq - 1)
average_volume = total_volume/total_weight

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

Iq = PqSq
Iq_beta = PqSq_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.

to:

1530834753223138

• pkienzle commented:

sasmodels/explore/beta/sasfit_compare.py in the beta_approx branch contains code for computing <F^2^> P S I=PS S_eff and I_beta = PS_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.

[RichardHeenan]

Trac update at 2018/08/22 16:05:01: richardh commented:

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

https://github.com/gregsuczewski/sasmodels

[RichardHeenan]

Trac update at 2018/09/28 16:30:30: richardh commented:

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-4volumevolfraction
http://trac.sasview.org/ticket/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-4volumevolfraction
http://trac.sasview.org/ticket/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.

[RichardHeenan]

Trac update at 2018/09/28 16:31:38: richardh changed attachment from "" to "sphere1.pdf"

[RichardHeenan]

Trac update at 2018/09/28 16:31:51: richardh changed attachment from "" to "sphere2.pdf"

[RichardHeenan]

Trac update at 2018/09/28 16:32:07: richardh changed attachment from "" to "sphere2a.pdf"

[RichardHeenan]

Trac update at 2018/09/28 16:32:44: richardh changed attachment from "" to "sphere3.pdf"

[RichardHeenan]

Trac update at 2018/09/28 16:32:56: richardh changed attachment from "" to "ellipsoid4..pdf"