source: sasmodels/explore/beta/sasfit_compare.py @ cddfef6

ticket_1156ticket_822_more_unit_tests
Last change on this file since cddfef6 was a34b811, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code
  • Property mode set to 100644
File size: 29.2 KB
Line 
1from __future__ import division, print_function
2# Make sasmodels available on the path
3import sys, os
4BETA_DIR = os.path.dirname(os.path.realpath(__file__))
5SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR))
6sys.path.insert(0, SASMODELS_DIR)
7
8from collections import namedtuple
9
10from matplotlib import pyplot as plt
11import numpy as np
12from numpy import pi, sin, cos, sqrt, fabs
13from numpy.polynomial.legendre import leggauss
14from scipy.special import j1 as J1
15from numpy import inf
16from scipy.special import gammaln  # type: ignore
17
18Theory = namedtuple('Theory', 'Q F1 F2 P S I Seff Ibeta')
19Theory.__new__.__defaults__ = (None,) * len(Theory._fields)
20
21#Used to calculate F(q) for the cylinder, sphere, ellipsoid models
22def sas_sinx_x(x):
23    with np.errstate(all='ignore'):
24        retvalue = sin(x)/x
25    retvalue[x == 0.] = 1.
26    return retvalue
27
28def sas_2J1x_x(x):
29    with np.errstate(all='ignore'):
30        retvalue = 2*J1(x)/x
31    retvalue[x == 0] = 1.
32    return retvalue
33
34def sas_3j1x_x(x):
35    """return 3*j1(x)/x"""
36    retvalue = np.empty_like(x)
37    with np.errstate(all='ignore'):
38        # GSL bessel_j1 taylor expansion
39        index = (x < 0.25)
40        y = x[index]**2
41        c1 = -1.0/10.0
42        c2 = +1.0/280.0
43        c3 = -1.0/15120.0
44        c4 = +1.0/1330560.0
45        c5 = -1.0/172972800.0
46        retvalue[index] = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))))
47        index = ~index
48        y = x[index]
49        retvalue[index] = 3*(sin(y) - y*cos(y))/y**3
50    retvalue[x == 0.] = 1.
51    return retvalue
52
53#Used to cross check my models with sasview models
54def build_model(model_name, q, **pars):
55    from sasmodels.core import load_model_info, build_model as build_sasmodel
56    from sasmodels.data import empty_data1D
57    from sasmodels.direct_model import DirectModel
58    model_info = load_model_info(model_name)
59    model = build_sasmodel(model_info, dtype='double!')
60    data = empty_data1D(q)
61    calculator = DirectModel(data, model,cutoff=0)
62    calculator.pars = pars.copy()
63    calculator.pars.setdefault('background', 0)
64    return calculator
65
66#gives the hardsphere structure factor that sasview uses
67def _hardsphere_simple(q, radius_effective, volfraction):
68    CUTOFFHS = 0.05
69    if fabs(radius_effective) < 1.E-12:
70        HARDSPH = 1.0
71        return HARDSPH
72    X = 1.0/(1.0 -volfraction)
73    D = X*X
74    A = (1.+2.*volfraction)*D
75    A *= A
76    X = fabs(q*radius_effective*2.0)
77    if X < 5.E-06:
78        HARDSPH = 1./A
79        return HARDSPH
80    X2 = X*X
81    B = (1.0 +0.5*volfraction)*D
82    B *= B
83    B *= -6.*volfraction
84    G = 0.5*volfraction*A
85    if X < CUTOFFHS:
86        FF = 8.0*A +6.0*B + 4.0*G + ( -0.8*A -B/1.5 -0.5*G +(A/35. +0.0125*B +0.02*G)*X2)*X2
87        HARDSPH = 1./(1. + volfraction*FF )
88        return HARDSPH
89    X4 = X2*X2
90    S, C = sin(X), cos(X)
91    FF =  ((G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S -(X2-2.)*C -2.) )/X + A*(S-X*C))/X
92    HARDSPH = 1./(1. + 24.*volfraction*FF/X2)
93    return HARDSPH
94
95def hardsphere_simple(q, radius_effective, volfraction):
96    SQ = [_hardsphere_simple(qk, radius_effective, volfraction) for qk in q]
97    return np.array(SQ)
98
99#Used in gaussian quadrature for polydispersity
100#returns values and the probability of those values based on gaussian distribution
101N_GAUSS = 35
102NSIGMA_GAUSS = 3
103def gaussian_distribution(center, sigma, lb, ub):
104    #3 standard deviations covers approx. 99.7%
105    if sigma != 0:
106        nsigmas = NSIGMA_GAUSS
107        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_GAUSS)
108        x = x[(x >= lb) & (x <= ub)]
109        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
110        return x, px
111    else:
112        return np.array([center]), np.array([1])
113
114N_SCHULZ = 80
115NSIGMA_SCHULZ = 8
116def schulz_distribution(center, sigma, lb, ub):
117    if sigma != 0:
118        nsigmas = NSIGMA_SCHULZ
119        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_SCHULZ)
120        x = x[(x >= lb) & (x <= ub)]
121        R = x/center
122        z = (center/sigma)**2
123        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
124        px = np.exp(arg)
125        return x, px
126    else:
127        return np.array([center]), np.array([1])
128
129#returns the effective radius used in sasview
130def ER_ellipsoid(radius_polar, radius_equatorial):
131    ee = np.empty_like(radius_polar)
132    if radius_polar > radius_equatorial:
133        ee = (radius_polar**2 - radius_equatorial**2)/radius_polar**2
134    elif radius_polar < radius_equatorial:
135        ee = (radius_equatorial**2 - radius_polar**2) / radius_equatorial**2
136    else:
137        ee = 2*radius_polar
138    if radius_polar * radius_equatorial != 0:
139        bd = 1.0 - ee
140        e1 = np.sqrt(ee)
141        b1 = 1.0 + np.arcsin(e1) / (e1*np.sqrt(bd))
142        bL = (1.0 + e1) / (1.0 - e1)
143        b2 = 1.0 + bd / 2 / e1 * np.log(bL)
144        delta = 0.75 * b1 * b2
145    ddd = np.zeros_like(radius_polar)
146    ddd = 2.0*(delta + 1.0)*radius_polar*radius_equatorial**2
147    return 0.5*ddd**(1.0 / 3.0)
148
149def ellipsoid_volume(radius_polar, radius_equatorial):
150    volume = (4./3.)*pi*radius_polar*radius_equatorial**2
151    return volume
152
153# F1 is F(q)
154# F2 is F(g)^2
155#IQM is I(q) with monodispersity
156#IQSM is I(q) with structure factor S(q) and monodispersity
157#IQBM is I(q) with Beta Approximation and monodispersity
158#SQ is monodisperse approach for structure factor
159#SQ_EFF is the effective structure factor from beta approx
160def ellipsoid_theta(q, radius_polar, radius_equatorial, sld, sld_solvent,
161                    volfraction=0, radius_effective=None):
162    #creates values z and corresponding probabilities w from legendre-gauss quadrature
163    volume = ellipsoid_volume(radius_polar, radius_equatorial)
164    z, w = leggauss(76)
165    F1 = np.zeros_like(q)
166    F2 = np.zeros_like(q)
167    #use a u subsition(u=cos) and then u=(z+1)/2 to change integration from
168    #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use
169    #legendre-gauss quadrature
170    for k, qk in enumerate(q):
171        r = sqrt(radius_equatorial**2*(1-((z+1)/2)**2)+radius_polar**2*((z+1)/2)**2)
172        form = (sld-sld_solvent)*volume*sas_3j1x_x(qk*r)
173        F2[k] = np.sum(w*form**2)
174        F1[k] = np.sum(w*form)
175    #the 1/2 comes from the change of variables mentioned above
176    F2 = F2/2.0
177    F1 = F1/2.0
178    if radius_effective is None:
179        radius_effective = ER_ellipsoid(radius_polar,radius_equatorial)
180    SQ = hardsphere_simple(q, radius_effective, volfraction)
181    SQ_EFF = 1 + F1**2/F2*(SQ - 1)
182    IQM = 1e-4*F2/volume
183    IQSM = volfraction*IQM*SQ
184    IQBM = volfraction*IQM*SQ_EFF
185    return Theory(Q=q, F1=F1, F2=F2, P=IQM, S=SQ, I=IQSM, Seff=SQ_EFF, Ibeta=IQBM)
186
187#IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc.
188#IQBD HAS NOT BEEN CROSS CHECKED AT ALL
189def ellipsoid_pe(q, radius_polar, radius_equatorial, sld, sld_solvent,
190                 radius_polar_pd=0.1, radius_equatorial_pd=0.1,
191                 radius_polar_pd_type='gaussian',
192                 radius_equatorial_pd_type='gaussian',
193                 volfraction=0, radius_effective=None,
194                 background=0, scale=1,
195                 norm='sasview'):
196    if norm not in ['sasview', 'sasfit', 'yun']:
197        raise TypeError("unknown norm "+norm)
198    if radius_polar_pd_type == 'gaussian':
199        Rp_val, Rp_prob = gaussian_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf)
200    elif radius_polar_pd_type == 'schulz':
201        Rp_val, Rp_prob = schulz_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf)
202    if radius_equatorial_pd_type == 'gaussian':
203        Re_val, Re_prob = gaussian_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf)
204    elif radius_equatorial_pd_type == 'schulz':
205        Re_val, Re_prob = schulz_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf)
206    total_weight = total_volume = 0
207    radius_eff = 0
208    F1, F2 = np.zeros_like(q), np.zeros_like(q)
209    for k, Rpk in enumerate(Rp_val):
210        print("ellipsoid cycle", k, "of", len(Rp_val))
211        for i, Rei in enumerate(Re_val):
212            theory = ellipsoid_theta(q, Rpk, Rei, sld, sld_solvent)
213            volume = ellipsoid_volume(Rpk, Rei)
214            weight = Rp_prob[k]*Re_prob[i]
215            total_weight += weight
216            total_volume += weight*volume
217            F1 += theory.F1*weight
218            F2 += theory.F2*weight
219            radius_eff += weight*ER_ellipsoid(Rpk, Rei)
220    F1 /= total_weight
221    F2 /= total_weight
222    average_volume = total_volume/total_weight
223    if radius_effective is None:
224        radius_effective = radius_eff/total_weight
225    if norm == 'sasfit':
226        IQD = F2
227    elif norm == 'sasview':
228        # Note: internally, sasview uses F2/total_volume because:
229        #   average_volume = total_volume/total_weight
230        #   F2/total_weight / average_volume
231        #     = F2/total_weight / total_volume/total_weight
232        #     = F2/total_volume
233        IQD = F2/average_volume*1e-4*volfraction
234        F1 *= 1e-2  # Yun is using sld in 1/A^2, not 1e-6/A^2
235        F2 *= 1e-4
236    elif norm == 'yun':
237        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
238        F2 *= 1e-12
239        IQD = F2/average_volume*1e8*volfraction
240    SQ = hardsphere_simple(q, radius_effective, volfraction)
241    beta = F1**2/F2
242    SQ_EFF = 1 + beta*(SQ - 1)
243    IQSD = IQD*SQ
244    IQBD = IQD*SQ_EFF
245    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
246
247#polydispersity for sphere
248def sphere_r(q,radius,sld,sld_solvent,
249             radius_pd=0.1, radius_pd_type='gaussian',
250             volfraction=0, radius_effective=None,
251             background=0, scale=1,
252             norm='sasview'):
253    if norm not in ['sasview', 'sasfit', 'yun']:
254        raise TypeError("unknown norm "+norm)
255    if radius_pd_type == 'gaussian':
256        radius_val, radius_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf)
257    elif radius_pd_type == 'schulz':
258        radius_val, radius_prob = schulz_distribution(radius, radius_pd*radius, 0, inf)
259    total_weight = total_volume = 0
260    F1 = np.zeros_like(q)
261    F2 = np.zeros_like(q)
262    for k, rk in enumerate(radius_val):
263        volume = 4./3.*pi*rk**3
264        total_weight += radius_prob[k]
265        total_volume += radius_prob[k]*volume
266        form = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk)
267        F2 += radius_prob[k]*form**2
268        F1 += radius_prob[k]*form
269    F1 /= total_weight
270    F2 /= total_weight
271    average_volume = total_volume/total_weight
272
273    if radius_effective is None:
274        radius_effective = radius
275    average_volume = total_volume/total_weight
276    if norm == 'sasfit':
277        IQD = F2
278    elif norm == 'sasview':
279        IQD = F2/average_volume*1e-4*volfraction
280    elif norm == 'yun':
281        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
282        F2 *= 1e-12
283        IQD = F2/average_volume*1e8*volfraction
284    SQ = hardsphere_simple(q, radius_effective, volfraction)
285    beta = F1**2/F2
286    SQ_EFF = 1 + beta*(SQ - 1)
287    IQSD = IQD*SQ
288    IQBD = IQD*SQ_EFF
289    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
290
291###############################################################################
292###############################################################################
293###############################################################################
294##################                                           ##################
295##################                   TESTS                   ##################
296##################                                           ##################
297###############################################################################
298###############################################################################
299###############################################################################
300
301def popn(d, keys):
302    """
303    Splits a dict into two, with any key of *d* which is in *keys* removed
304    from *d* and added to *b*. Returns *b*.
305    """
306    b = {}
307    for k in keys:
308        try:
309            b[k] = d.pop(k)
310        except KeyError:
311            pass
312    return b
313
314def sasmodels_theory(q, Pname, **pars):
315    """
316    Call sasmodels to compute the model with and without a hard sphere
317    structure factor.
318    """
319    #mono_pars = {k: (0 if k.endswith('_pd') else v) for k, v in pars.items()}
320    Ppars = pars.copy()
321    Spars = popn(Ppars, ['radius_effective', 'volfraction'])
322    Ipars = pars.copy()
323
324    # Autofill npts and nsigmas for the given pd type
325    for k, v in pars.items():
326        if k.endswith("_pd_type"):
327            if v == "gaussian":
328                n, nsigmas = N_GAUSS, NSIGMA_GAUSS
329            elif v == "schulz":
330                n, nsigmas = N_SCHULZ, NSIGMA_SCHULZ
331            Ppars.setdefault(k.replace("_pd_type", "_pd_n"), n)
332            Ppars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas)
333            Ipars.setdefault(k.replace("_pd_type", "_pd_n"), n)
334            Ipars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas)
335
336    #Ppars['scale'] = Spars.get('volfraction', 1)
337    P = build_model(Pname, q)
338    S = build_model("hardsphere", q)
339    I = build_model(Pname+"@hardsphere", q)
340    Pq = P(**Ppars)*pars.get('volfraction', 1)
341    Sq = S(**Spars)
342    Iq = I(**Ipars)
343    #Iq = Pq*Sq*pars.get('volfraction', 1)
344    #Sq = Iq/Pq
345    #Iq = None#= Sq = None
346    r = dict(I._kernel.results())
347    return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r["S_eff(Q)"], Ibeta=Iq)
348
349def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'):
350    """
351    Plot fields in common between target and actual, along with relative error.
352    """
353    available = [s for s in fields.split()
354                 if getattr(target, s) is not None and getattr(actual, s) is not None]
355    rows = len(available)
356    for row, field in enumerate(available):
357        Q = target.Q
358        I1, I2 = getattr(target, field), getattr(actual, field)
359        plt.subplot(rows, 2, 2*row+1)
360        plt.loglog(Q, abs(I1), label="target "+field)
361        plt.loglog(Q, abs(I2), label="value "+field)
362        #plt.legend(loc="upper left", bbox_to_anchor=(1,1))
363        plt.legend(loc='lower left')
364        plt.subplot(rows, 2, 2*row+2)
365        plt.semilogx(Q, I2/I1 - 1, label="relative error")
366        #plt.semilogx(Q, I1/I2 - 1, label="relative error")
367    plt.tight_layout()
368    plt.suptitle(title)
369    plt.show()
370
371def data_file(name):
372    return os.path.join(BETA_DIR, 'data_files', name)
373
374def load_sasfit(path):
375    data = np.loadtxt(path, dtype=str, delimiter=';').T
376    data = np.vstack((map(float, v) for v in data[0:2]))
377    return data
378
379COMPARISON = {}  # Type: Dict[(str,str,str)] -> Callable[(), None]
380
381def compare_sasview_sphere(pd_type='schulz'):
382    q = np.logspace(-5, 0, 250)
383    model = 'sphere'
384    pars = dict(
385        radius=20, sld=4, sld_solvent=1,
386        background=0,
387        radius_pd=.1, radius_pd_type=pd_type,
388        volfraction=0.15,
389        #radius_effective=12.59921049894873,  # equivalent average sphere radius
390        )
391    target = sasmodels_theory(q, model, **pars)
392    actual = sphere_r(q, norm='sasview', **pars)
393    title = " ".join(("sasmodels", model, pd_type))
394    compare(title, target, actual)
395COMPARISON[('sasview', 'sphere', 'gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian')
396COMPARISON[('sasview', 'sphere', 'schulz')] = lambda: compare_sasview_sphere(pd_type='schulz')
397
398def compare_sasview_ellipsoid(pd_type='gaussian'):
399    q = np.logspace(-5, 0, 50)
400    model = 'ellipsoid'
401    pars = dict(
402        radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1,
403        background=0,
404        radius_polar_pd=0.1, radius_polar_pd_type=pd_type,
405        radius_equatorial_pd=0.1, radius_equatorial_pd_type=pd_type,
406        volfraction=0.15,
407        radius_effective=270.7543927018,
408        #radius_effective=12.59921049894873,
409        )
410    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars)
411    actual = ellipsoid_pe(q, norm='sasview', **pars)
412    title = " ".join(("sasmodels", model, pd_type))
413    compare(title, target, actual)
414COMPARISON[('sasview', 'ellipsoid', 'gaussian')] = lambda: compare_sasview_ellipsoid(pd_type='gaussian')
415COMPARISON[('sasview', 'ellipsoid', 'schulz')] = lambda: compare_sasview_ellipsoid(pd_type='schulz')
416
417def compare_yun_ellipsoid_mono():
418    pars = {
419        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
420        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
421        'sld': 2, 'sld_solvent': 1,
422        'volfraction': 0.15,
423        # Yun uses radius for same volume sphere for effective radius
424        # whereas sasview uses the average curvature.
425        'radius_effective': 12.59921049894873,
426    }
427
428    data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T
429    Q = data[0]
430    F1 = data[1]
431    P = data[3]*pars['volfraction']
432    S = data[5]
433    Seff = data[6]
434    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
435    actual = ellipsoid_pe(Q, norm='yun', **pars)
436    title = " ".join(("yun", "ellipsoid", "no pd"))
437    #compare(title, target, actual, fields="P S I Seff Ibeta")
438    compare(title, target, actual)
439COMPARISON[('yun', 'ellipsoid', 'gaussian')] = compare_yun_ellipsoid_mono
440COMPARISON[('yun', 'ellipsoid', 'schulz')] = compare_yun_ellipsoid_mono
441
442def compare_yun_sphere_gauss():
443    # Note: yun uses gauss limits from R0/10 to R0 + 5 sigma steps sigma/100
444    # With pd = 0.1, that's 14 sigma and 1400 points.
445    pars = {
446        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
447        'sld': 6, 'sld_solvent': 0,
448        'volfraction': 0.1,
449    }
450
451    data = np.loadtxt(data_file('testPolydisperseGaussianSphere.dat'), skiprows=2).T
452    Q = data[0]
453    F1 = data[1]
454    F2 = data[2]
455    P = data[3]
456    S = data[5]
457    Seff = data[6]
458    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
459    actual = sphere_r(Q, norm='yun', **pars)
460    title = " ".join(("yun", "sphere", "10% dispersion 10% Vf"))
461    compare(title, target, actual)
462    data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'), skiprows=2).T
463    pars.update(radius_pd=0.15)
464    Q = data[0]
465    F1 = data[1]
466    F2 = data[2]
467    P = data[3]
468    S = data[5]
469    Seff = data[6]
470    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
471    actual = sphere_r(Q, norm='yun', **pars)
472    title = " ".join(("yun", "sphere", "15% dispersion 10% Vf"))
473    compare(title, target, actual)
474COMPARISON[('yun', 'sphere', 'gaussian')] = compare_yun_sphere_gauss
475
476
477def compare_sasfit_sphere_gauss():
478    #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3
479    pars = {
480        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
481        'sld': 4, 'sld_solvent': 1,
482        'volfraction': 0.3,
483    }
484
485    Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt'))
486    Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt'))
487    Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt'))
488    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt'))
489    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt'))
490    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
491    actual = sphere_r(Q, norm="sasfit", **pars)
492    title = " ".join(("sasfit", "sphere", "pd=10% gaussian"))
493    compare(title, target, actual)
494    #compare(title, target, actual, fields="P")
495COMPARISON[('sasfit', 'sphere', 'gaussian')] = compare_sasfit_sphere_gauss
496
497def compare_sasfit_sphere_schulz():
498    #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1
499    #We have scaled the output from sasfit by 1e-4*volume*volfraction
500    #0.10050378152592121
501    pars = {
502        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz',
503        'sld': 4, 'sld_solvent': 1,
504        'volfraction': 0.3,
505    }
506
507    Q, IQ = load_sasfit(data_file('sasfit_sphere_schulz_IQD.txt'))
508    Q, IQSD = load_sasfit(data_file('sasfit_sphere_schulz_IQSD.txt'))
509    Q, IQBD = load_sasfit(data_file('sasfit_sphere_schulz_IQBD.txt'))
510    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
511    actual = sphere_r(Q, norm="sasfit", **pars)
512    title = " ".join(("sasfit", "sphere", "pd=10% schulz"))
513    compare(title, target, actual)
514COMPARISON[('sasfit', 'sphere', 'schulz')] = compare_sasfit_sphere_schulz
515
516def compare_sasfit_ellipsoid_schulz():
517    #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1
518    #Effective radius =13.1353356684
519    #We have scaled the output from sasfit by 1e-4*volume*volfraction
520    #0.10050378152592121
521    pars = {
522        'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz',
523        'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz',
524        'sld': 4, 'sld_solvent': 1,
525        'volfraction': 0.3, 'radius_effective': 13.1353356684,
526    }
527
528    Q, IQ = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQD.txt'))
529    Q, IQSD = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQSD.txt'))
530    Q, IQBD = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQBD.txt'))
531    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
532    actual = ellipsoid_pe(Q, norm="sasfit", **pars)
533    title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz"))
534    compare(title, target, actual)
535COMPARISON[('sasfit', 'ellipsoid', 'schulz')] = compare_sasfit_ellipsoid_schulz
536
537
538def compare_sasfit_ellipsoid_gaussian():
539    pars = {
540        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
541        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
542        'sld': 4, 'sld_solvent': 1,
543        'volfraction': 0, 'radius_effective': None,
544    }
545
546    #Rp=20,Re=10,eta_core=4,eta_solv=1
547    Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt'))
548    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None)
549    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
550    target = Theory(Q=Q, P=PQ0)
551    compare("sasfit ellipsoid no poly", target, actual); plt.show()
552
553    #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
554    Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt'))
555    pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None)
556    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
557    target = Theory(Q=Q, P=PQ_Rp10)
558    compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show()
559    #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
560    Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt'))
561    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None)
562    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
563    target = Theory(Q=Q, P=PQ_Re10)
564    compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show()
565    #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
566    Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt'))
567    pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None)
568    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
569    target = Theory(Q=Q, P=PQ_Rp30)
570    compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show()
571    #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
572    Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt'))
573    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None)
574    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
575    target = Theory(Q=Q, P=PQ_Re30)
576    compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show()
577    #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
578    Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt'))
579    pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None)
580    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
581    target = Theory(Q=Q, P=PQ_Rp60)
582    compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show()
583    #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
584    Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt'))
585    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None)
586    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
587    target = Theory(Q=Q, P=PQ_Re60)
588    compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show()
589
590    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15
591    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt'))
592    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt'))
593    pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
594    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
595    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
596    compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show()
597    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15
598    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt'))
599    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt'))
600    pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
601    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
602    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
603    compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show()
604    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15
605    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt'))
606    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt'))
607    pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
608    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
609    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
610    compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show()
611
612    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3
613    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt'))
614    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt'))
615    pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
616    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
617    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
618    compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show()
619    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3
620    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt'))
621    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt'))
622    pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
623    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
624    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
625    compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show()
626    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3
627    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt'))
628    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt'))
629    pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
630    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
631    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
632    compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show()
633
634    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6
635    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt'))
636    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt'))
637    pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
638    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
639    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
640    compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show()
641    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6
642    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt'))
643    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt'))
644    pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
645    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
646    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
647    compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show()
648    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6
649    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt'))
650    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt'))
651    pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
652    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
653    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
654    compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show()
655COMPARISON[('sasfit', 'ellipsoid', 'gaussian')] = compare_sasfit_ellipsoid_gaussian
656
657def main():
658    key = tuple(sys.argv[1:])
659    if key not in COMPARISON:
660        print("usage: sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid] [gaussian|schulz]")
661        return
662    comparison = COMPARISON[key]
663    comparison()
664
665if __name__ == "__main__":
666    main()
Note: See TracBrowser for help on using the repository browser.