Changeset 01c8d9e in sasmodels


Ignore:
Timestamp:
Aug 7, 2018 2:32:18 PM (4 months ago)
Author:
Suczewski <ges3@…>
Branches:
F1F2models_grethe, beta_approx, py3, ticket-1015-gpu-mem-error, ticket-1157, ticket-608-user-defined-weights, ticket_1156
Children:
7420664
Parents:
707cbdb
Message:

beta approximation, first pass

Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed
  • explore/beta/sasfit_compare.py

    r707cbdb r01c8d9e  
    33import sys,os 
    44BETA_DIR = os.path.dirname(os.path.realpath(__file__)) 
    5 SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) 
     5#SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) 
     6SASMODELS_DIR = r"C:\Source\sasmodels" 
    67sys.path.insert(0, SASMODELS_DIR) 
    78import os 
     9from collections import namedtuple 
     10 
    811from matplotlib import pyplot as plt 
    912import numpy as np 
     
    1417from scipy.special import gammaln  # type: ignore 
    1518 
    16 # THE FOLLOWING 7 BOOLEANS TAILOR WHICH GRAPHS ARE PRINTED WHEN THE PROGRAM RUNS 
    17 #RICHARD, YOU WILL WANT SASVIEW, SASFIT, SCHULZ, ELLIPSOID, AND SPHERE ALL TRUE. 
    18 ELLIPSOID = False 
    19 SPHERE = True 
    20  
    21 GAUSSIAN = True 
    22 SCHULZ = False 
    23  
    24 SASVIEW=True 
    25 SASFIT=True 
    26 YUN = True 
    27  
    28 def data_file(name): 
    29     return os.path.join(BETA_DIR + '\\data_files', name) 
     19Theory = namedtuple('Theory', 'Q F1 F2 P S I Seff Ibeta') 
     20Theory.__new__.__defaults__ = (None,) * len(Theory._fields) 
    3021 
    3122#Used to calculate F(q) for the cylinder, sphere, ellipsoid models 
     
    4738    with np.errstate(all='ignore'): 
    4839        # GSL bessel_j1 taylor expansion 
    49         index = (x < 0.25)   
     40        index = (x < 0.25) 
    5041        y = x[index]**2 
    5142        c1 = -1.0/10.0 
     
    7162    calculator = DirectModel(data, model,cutoff=0) 
    7263    calculator.pars = pars.copy() 
    73     calculator.pars.setdefault('background', 1e-3) 
     64    calculator.pars.setdefault('background', 0) 
    7465    return calculator 
    7566 
    7667#gives the hardsphere structure factor that sasview uses 
    77 def hardsphere_simple(q, radius_effective, volfraction):  
    78     CUTOFFHS=0.05  
     68def _hardsphere_simple(q, radius_effective, volfraction): 
     69    CUTOFFHS=0.05 
    7970    if fabs(radius_effective) < 1.E-12: 
    8071        HARDSPH=1.0 
     
    9687        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 
    9788        HARDSPH= 1./(1. + volfraction*FF ) 
    98         return HARDSPH    
     89        return HARDSPH 
    9990    X4=X2*X2 
    10091    S, C = sin(X), cos(X) 
    101     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 ; 
    102     HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); 
     92    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 
     93    HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ) 
    10394    return HARDSPH 
     95 
     96def hardsphere_simple(q, radius_effective, volfraction): 
     97    SQ = [_hardsphere_simple(qk, radius_effective, volfraction) for qk in q] 
     98    return np.array(SQ) 
    10499 
    105100#Used in gaussian quadrature for polydispersity 
    106101#returns values and the probability of those values based on gaussian distribution 
    107 def gaussian_distribution(center, sigma,lb,ub): 
    108     #3 standard deviations covers approx. 99.7%  
     102N_GAUSS = 35 
     103NSIGMA_GAUSS = 3 
     104def gaussian_distribution(center, sigma, lb, ub): 
     105    #3 standard deviations covers approx. 99.7% 
    109106    if sigma != 0: 
    110         nsigmas=3 
    111         x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=35) 
     107        nsigmas = NSIGMA_GAUSS 
     108        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_GAUSS) 
    112109        x= x[(x >= lb) & (x <= ub)] 
    113110        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma)) 
     
    116113        return np.array([center]), np.array([1]) 
    117114 
     115N_SCHULZ = 80 
     116NSIGMA_SCHULZ = 8 
    118117def schulz_distribution(center, sigma, lb, ub): 
    119118    if sigma != 0: 
    120         nsigmas=8 
    121         x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=80) 
     119        nsigmas = NSIGMA_SCHULZ 
     120        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_SCHULZ) 
    122121        x= x[(x >= lb) & (x <= ub)] 
    123122        R = x/center 
     
    160159#SQ is monodisperse approach for structure factor 
    161160#SQ_EFF is the effective structure factor from beta approx 
    162 def ellipsoid_Theta(q, radius_polar, radius_equatorial, sld, sld_solvent,volfraction=0,effective_radius=0): 
     161def ellipsoid_theta(q, radius_polar, radius_equatorial, sld, sld_solvent, 
     162                    volfraction=0, radius_effective=None): 
    163163    #creates values z and corresponding probabilities w from legendre-gauss quadrature 
     164    volume = ellipsoid_volume(radius_polar, radius_equatorial) 
    164165    z, w = leggauss(76) 
    165166    F1 = np.zeros_like(q) 
    166167    F2 = np.zeros_like(q) 
    167168    #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    #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use 
    169170    #legendre-gauss quadrature 
    170171    for k, qk in enumerate(q): 
    171172        r = sqrt(radius_equatorial**2*(1-((z+1)/2)**2)+radius_polar**2*((z+1)/2)**2) 
    172         F2i = ((sld-sld_solvent)*ellipsoid_volume(radius_polar,radius_equatorial)*sas_3j1x_x(qk*r))**2 
    173         F2[k] = np.sum(w*F2i) 
    174         F1i = (sld-sld_solvent)*ellipsoid_volume(radius_polar,radius_equatorial)*sas_3j1x_x(qk*r) 
    175         F1[k] = np.sum(w*F1i) 
     173        form = (sld-sld_solvent)*volume*sas_3j1x_x(qk*r) 
     174        F2[k] = np.sum(w*form**2) 
     175        F1[k] = np.sum(w*form) 
    176176    #the 1/2 comes from the change of variables mentioned above 
    177177    F2 = F2/2.0 
    178178    F1 = F1/2.0 
    179     if effective_radius == 0: 
    180         effective_radius = ER_ellipsoid(radius_polar,radius_equatorial) 
    181     else: 
    182         effective_radius = effective_radius 
    183     SQ = np.array([hardsphere_simple(qk, effective_radius, volfraction) for qk in q])    
    184     SQ_EFF = 1 + F1**2/F2*(SQ - 1)  
    185     IQM = 1e-4/ellipsoid_volume(radius_polar,radius_equatorial)*F2 
     179    if radius_effective is None: 
     180        radius_effective = ER_ellipsoid(radius_polar,radius_equatorial) 
     181    SQ = hardsphere_simple(q, radius_effective, volfraction) 
     182    SQ_EFF = 1 + F1**2/F2*(SQ - 1) 
     183    IQM = 1e-4*F2/volume 
    186184    IQSM = volfraction*IQM*SQ 
    187185    IQBM = volfraction*IQM*SQ_EFF 
    188     return F1, F2, IQM, IQSM, IQBM, SQ, SQ_EFF 
    189  
    190 #IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc.  
     186    return Theory(Q=q, F1=F1, F2=F2, P=IQM, S=SQ, I=IQSM, Seff=SQ_EFF, Ibeta=IQBM) 
     187 
     188#IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc. 
    191189#IQBD HAS NOT BEEN CROSS CHECKED AT ALL 
    192 def ellipsoid_pe(q, Rp, Re, sld, sld_solvent, volfraction=0,effective_radius=0,radius_polar_pd=0.1,radius_equatorial_pd=0.1,distribution='gaussian'): 
    193     if distribution == 'gaussian': 
    194         Rp_val, Rp_prob = gaussian_distribution(Rp, radius_polar_pd*Rp, 0, inf) 
    195         Re_val, Re_prob = gaussian_distribution(Re, radius_equatorial_pd*Re, 0, inf) 
    196     elif distribution == 'schulz': 
    197         Rp_val, Rp_prob = schulz_distribution(Rp, radius_polar_pd*Rp, 0, inf) 
    198         Re_val, Re_prob = schulz_distribution(Re, radius_equatorial_pd*Re, 0, inf) 
    199     Normalization = 0 
    200     total_weight = 0 
    201     PQ = np.zeros_like(q) 
    202     F12,F21 = np.zeros_like(q), np.zeros_like(q) 
     190def ellipsoid_pe(q, radius_polar, radius_equatorial, sld, sld_solvent, 
     191                 radius_polar_pd=0.1, radius_equatorial_pd=0.1, 
     192                 radius_polar_pd_type='gaussian', 
     193                 radius_equatorial_pd_type='gaussian', 
     194                 volfraction=0, radius_effective=None, 
     195                 background=0, scale=1, 
     196                 norm='sasview'): 
     197    if norm not in ['sasview', 'sasfit', 'yun']: 
     198        raise TypeError("unknown norm "+norm) 
     199    if radius_polar_pd_type == 'gaussian': 
     200        Rp_val, Rp_prob = gaussian_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf) 
     201    elif radius_polar_pd_type == 'schulz': 
     202        Rp_val, Rp_prob = schulz_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf) 
     203    if radius_equatorial_pd_type == 'gaussian': 
     204        Re_val, Re_prob = gaussian_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf) 
     205    elif radius_equatorial_pd_type == 'schulz': 
     206        Re_val, Re_prob = schulz_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf) 
     207    total_weight = total_volume = 0 
    203208    radius_eff = 0 
     209    F1, F2 = np.zeros_like(q), np.zeros_like(q) 
    204210    for k, Rpk in enumerate(Rp_val): 
    205211        for i, Rei in enumerate(Re_val): 
    206             F1i, F2i, PQi, IQSM, IQBM, SQ, SQ_EFF = ellipsoid_Theta(q,Rpk,Rei,sld,sld_solvent) 
    207             radius_eff += Rp_prob[k]*Re_prob[i]*ER_ellipsoid(Rpk,Rei) 
    208             total_weight +=  Rp_prob[k]*Re_prob[i] 
    209             Normalization += Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei) 
    210             PQ += PQi*Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei) 
    211             F12 += F1i*Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei) 
    212             F21 += F2i*Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei) 
    213     F12 = (F12/Normalization)**2 
    214     F21 = F21/Normalization 
    215     if effective_radius == 0: 
    216         effective_radius = radius_eff/total_weight 
    217     else: 
    218         effective_radius = effective_radius 
    219     SQ = np.array([hardsphere_simple(qk, effective_radius, volfraction) for qk in q])    
    220     SQ_EFF = 1 + F12/F21*(SQ - 1)      
    221     IQD = PQ/Normalization 
    222     IQSD = volfraction*IQD*SQ 
    223     IQBD = volfraction*IQD*SQ_EFF 
    224     return IQD, IQSD, IQBD, SQ, SQ_EFF 
     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    print("this is the effective radius for pure python",radius_effective) 
     226    if norm == 'sasfit': 
     227        IQD = F2 
     228    elif norm == 'sasview': 
     229        # Note: internally, sasview uses F2/total_volume because: 
     230        #   average_volume = total_volume/total_weight 
     231        #   F2/total_weight / average_volume 
     232        #     = F2/total_weight / total_volume/total_weight 
     233        #     = F2/total_volume 
     234         
     235        IQD = F2/average_volume*1e-4*volfraction 
     236        F1 *= 1e-2  # Yun is using sld in 1/A^2, not 1e-6/A^2 
     237        F2 *= 1e-4 
     238    elif norm == 'yun': 
     239        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2 
     240        F2 *= 1e-12 
     241        IQD = F2/average_volume*1e8*volfraction 
     242    SQ = hardsphere_simple(q, radius_effective, volfraction) 
     243    beta = F1**2/F2 
     244    SQ_EFF = 1 + beta*(SQ - 1) 
     245    IQSD = IQD*SQ 
     246    IQBD = IQD*SQ_EFF 
     247    print("\n\n\n\n\n this is F1" + str(F1)) 
     248 
     249    print("\n\n\n\n\n this is F2" + str(F2)) 
     250    print("\n\n\n\n\n this is SQ" + str(SQ))    
     251    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) 
    225252 
    226253#polydispersity for sphere 
    227 def sphere_r(q,radius,sld,sld_solvent,volfraction=0,radius_pd=0.1,distribution='gaussian',norm_type='volfraction'): 
    228     if distribution == 'gaussian':     
     254def sphere_r(q,radius,sld,sld_solvent, 
     255             radius_pd=0.1, radius_pd_type='gaussian', 
     256             volfraction=0, radius_effective=None, 
     257             background=0, scale=1, 
     258             norm='sasview'): 
     259    if norm not in ['sasview', 'sasfit', 'yun']: 
     260        raise TypeError("unknown norm "+norm) 
     261    if radius_pd_type == 'gaussian': 
    229262        radius_val, radius_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf) 
    230     elif distribution == 'schulz': 
     263    elif radius_pd_type == 'schulz': 
    231264        radius_val, radius_prob = schulz_distribution(radius, radius_pd*radius, 0, inf) 
    232     Normalization=0 
     265    total_weight = total_volume = 0 
     266    F1 = np.zeros_like(q) 
    233267    F2 = np.zeros_like(q) 
    234     F1 = np.zeros_like(q) 
    235     if norm_type == 'numdensity': 
    236         for k, rk in enumerate(radius_val): 
    237             volume = 4./3.*pi*rk**3 
    238             Normalization += radius_prob[k]*volume 
    239             F2i = ((sld-sld_solvent)*volume*sas_3j1x_x(q*rk))**2/volume 
    240             F1i = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk)/volume 
    241              
    242             F2 += radius_prob[k]*volume*F2i 
    243             F1 += radius_prob[k]*volume*F1i 
    244         F21 = F2/Normalization 
    245         F12 = (F1/Normalization)**2 
    246         SQ = np.array([hardsphere_simple(qk, radius, volfraction) for qk in q]) 
    247         SQ_EFF = 1 + F12/F21*(SQ - 1)  
    248         IQD = 1e-4*F21 
    249         IQSD = volfraction*IQD*SQ 
    250         IQBD = volfraction*IQD*SQ_EFF 
    251     elif norm_type == 'volfraction': 
    252         for k, rk in enumerate(radius_val): 
    253             volume = 4./3.*pi*rk**3 
    254             Normalization += radius_prob[k] 
    255             F2i = ((sld-sld_solvent)*volume*sas_3j1x_x(q*rk))**2 
    256             F1i = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk) 
    257             F2 += radius_prob[k]*F2i 
    258             F1 += radius_prob[k]*F1i 
    259      
    260         F21 = F2/Normalization 
    261         F12 = (F1/Normalization)**2 
    262         SQ = np.array([hardsphere_simple(qk, radius, volfraction) for qk in q]) 
    263         SQ_EFF = 1 + F12/F21*(SQ - 1)  
    264         IQD = 1e-4/(4./3.*pi*radius**3)*F21 
    265         IQSD = volfraction*IQD*SQ 
    266         IQBD = volfraction*IQD*SQ_EFF 
    267     return IQD, IQSD, IQBD, SQ, SQ_EFF 
     268    for k, rk in enumerate(radius_val): 
     269        volume = 4./3.*pi*rk**3 
     270        total_weight += radius_prob[k] 
     271        total_volume += radius_prob[k]*volume 
     272        form = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk) 
     273        F2 += radius_prob[k]*form**2 
     274        F1 += radius_prob[k]*form 
     275    F1 /= total_weight 
     276    F2 /= total_weight 
     277    average_volume = total_volume/total_weight 
     278 
     279    if radius_effective is None: 
     280        radius_effective = radius 
     281    average_volume = total_volume/total_weight 
     282    if norm == 'sasfit': 
     283        IQD = F2 
     284    elif norm == 'sasview': 
     285        IQD = F2/average_volume*1e-4*volfraction 
     286    elif norm == 'yun': 
     287        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2 
     288        F2 *= 1e-12 
     289        IQD = F2/average_volume*1e8*volfraction 
     290    SQ = hardsphere_simple(q, radius_effective, volfraction) 
     291    beta = F1**2/F2 
     292    SQ_EFF = 1 + beta*(SQ - 1) 
     293    IQSD = IQD*SQ 
     294    IQBD = IQD*SQ_EFF 
     295    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) 
    268296 
    269297############################################################################### 
     
    276304############################################################################### 
    277305############################################################################### 
    278 # SASVIEW     
    279 if (ELLIPSOID == True) & (GAUSSIAN == True) & (SASVIEW == True): 
     306 
     307def popn(d, keys): 
     308    """ 
     309    Splits a dict into two, with any key of *d* which is in *keys* removed 
     310    from *d* and added to *b*. Returns *b*. 
     311    """ 
     312    b = {} 
     313    for k in keys: 
     314        try: 
     315            b[k] = d.pop(k) 
     316        except KeyError: 
     317            pass 
     318    return b 
     319 
     320def sasmodels_theory(q, Pname, **pars): 
     321    """ 
     322    Call sasmodels to compute the model with and without a hard sphere 
     323    structure factor. 
     324    """ 
     325    #mono_pars = {k: (0 if k.endswith('_pd') else v) for k, v in pars.items()} 
     326    Ppars = pars.copy() 
     327    Spars = popn(Ppars, ['radius_effective', 'volfraction']) 
     328    Ipars = pars.copy() 
     329 
     330    # Autofill npts and nsigmas for the given pd type 
     331    for k, v in pars.items(): 
     332        if k.endswith("_pd_type"): 
     333            if v == "gaussian": 
     334                n, nsigmas = N_GAUSS, NSIGMA_GAUSS 
     335            elif v == "schulz": 
     336                n, nsigmas = N_SCHULZ, NSIGMA_SCHULZ 
     337            Ppars.setdefault(k.replace("_pd_type", "_pd_n"), n) 
     338            Ppars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas) 
     339            Ipars.setdefault(k.replace("_pd_type", "_pd_n"), n) 
     340            Ipars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas) 
     341 
     342     
     343    #Ppars['scale'] = Spars.get('volfraction', 1) 
     344    P = build_model(Pname, q) 
     345    S = build_model("hardsphere", q) 
     346    I = build_model(Pname+"@hardsphere", q) 
     347    Pq = P(**Ppars)*pars.get('volfraction', 1) 
     348    Sq = S(**Spars) 
     349    Iq = I(**Ipars) 
     350    #Iq = Pq*Sq*pars.get('volfraction', 1) 
     351    #Sq = Iq/Pq 
     352    #Iq = None#= Sq = None 
     353    r=I._kernel.results 
     354    return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r[1], Ibeta=Iq) 
     355 
     356def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'): 
     357    """ 
     358    Plot fields in common between target and actual, along with relative error. 
     359    """ 
     360    available = [s for s in fields.split() 
     361                 if getattr(target, s) is not None and getattr(actual, s) is not None] 
     362    rows = len(available) 
     363    for row, field in enumerate(available): 
     364        Q = target.Q 
     365        I1, I2 = getattr(target, field), getattr(actual, field) 
     366        plt.subplot(rows, 2, 2*row+1) 
     367        plt.loglog(Q, abs(I1), label="target "+field) 
     368        plt.loglog(Q, abs(I2), label="value "+field) 
     369        #plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
     370        plt.legend(loc='lower left') 
     371        plt.subplot(rows, 2, 2*row+2) 
     372        plt.semilogx(Q, I2/I1 - 1, label="relative error") 
     373        #plt.semilogx(Q, I1/I2 - 1, label="relative error") 
     374    plt.tight_layout() 
     375    plt.suptitle(title) 
     376    plt.show() 
     377 
     378def data_file(name): 
     379    return os.path.join(BETA_DIR, 'data_files', name) 
     380 
     381def load_sasfit(path): 
     382    data = np.loadtxt(path, dtype=str, delimiter=';').T 
     383    data = np.vstack((map(float, v) for v in data[0:2])) 
     384    return data 
     385 
     386COMPARISON = {}  # Type: Dict[(str,str,str)] -> Callable[(), None] 
     387 
     388def compare_sasview_sphere(pd_type='schulz'): 
     389    q = np.logspace(-5, 0, 250) 
     390    model = 'sphere' 
     391    pars = dict( 
     392        radius=20,sld=4,sld_solvent=1, 
     393        background=0, 
     394        radius_pd=.1, radius_pd_type=pd_type, 
     395        volfraction=0.15, 
     396        #radius_effective=12.59921049894873,  # equivalent average sphere radius 
     397        ) 
     398    target = sasmodels_theory(q, model, **pars) 
     399    actual = sphere_r(q, norm='sasview', **pars) 
     400    title = " ".join(("sasmodels", model, pd_type)) 
     401    compare(title, target, actual) 
     402COMPARISON[('sasview','sphere','gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian') 
     403COMPARISON[('sasview','sphere','schulz')] = lambda: compare_sasview_sphere(pd_type='schulz') 
     404 
     405def compare_sasview_ellipsoid(pd_type='gaussian'): 
    280406    q = np.logspace(-5, 0, 50) 
    281     volfraction = 0.2 
    282     model = build_model("ellipsoid",q) 
    283     Sq = model(radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1,\ 
    284                background=0,radius_polar_pd=.1,radius_polar_pd_n=35,radius_equatorial_pd=.1,radius_equatorial_pd_n=35) 
    285     IQD, IQSD, IQBD, SQ, SQ_EFF = ellipsoid_pe(q,20,400,4,1,volfraction) 
    286     Sq2 = model(radius_polar=20,radius_equatorial=10,sld=2,sld_solvent=1,background=0) 
    287     IQM = ellipsoid_Theta(q,20,10,2,1)[2] 
    288     model2 = build_model("ellipsoid@hardsphere", q) 
    289     Sq3 = model2(radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1,background=0,radius_polar_pd=.1,\ 
    290                  radius_polar_pd_n=35,radius_equatorial_pd=.1,radius_equatorial_pd_n=35) 
    291     Sqp = model(radius_polar=20,radius_equatorial=10,sld=4,sld_solvent=1,\ 
    292                    background=0,radius_polar_pd=0.1,radius_polar_pd_n=35,radius_equatorial_pd=0,radius_equatorial_pd_n=35) 
    293     IQDp = ellipsoid_pe(q,20,10,4,1,radius_polar_pd=0.1,radius_equatorial_pd=0)[0]   
    294     Sqp2 = model(radius_polar=20,radius_equatorial=10,sld=4,sld_solvent=1,\ 
    295                    background=0,radius_polar_pd=0,radius_polar_pd_n=35,radius_equatorial_pd=0.1,radius_equatorial_pd_n=35) 
    296     IQDe = ellipsoid_pe(q,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.1)[0] 
    297     print('\n\n\n\n ELLIPSOID COMPARISON WITH SASVIEW WITHOUT V') 
    298     plt.subplot(323) 
    299     plt.loglog(q, Sq,'r') 
    300     plt.loglog(q,IQD,'b') 
    301     plt.title('IQD') 
    302     plt.subplot(324) 
    303     plt.semilogx(q,Sq/IQD-1) 
    304     plt.subplot(321) 
    305     plt.title('IQM') 
    306     plt.loglog(q, Sq2,'r') 
    307     plt.loglog(q,IQM,'b') 
    308     plt.subplot(322) 
    309     plt.semilogx(q,IQM/Sq2-1) 
    310     plt.subplot(325) 
    311     plt.title('IQSD') 
    312     plt.loglog(q, Sq3, '-b') 
    313     plt.loglog(q, IQSD, '-r') 
    314     plt.subplot(326) 
    315     plt.semilogx(q, Sq3/IQSD-1, 'b') 
    316     plt.tight_layout() 
    317     plt.show() 
    318     plt.subplot(221) 
    319     plt.loglog(q,IQDp,'r') 
    320     plt.loglog(q,Sqp,'b') 
    321     plt.title('IQD Polar') 
    322     plt.subplot(222) 
    323     plt.plot(q,IQDp/Sqp-1) 
    324     plt.subplot(223) 
    325     plt.loglog(q,IQDe,'r') 
    326     plt.loglog(q,Sqp2,'b') 
    327     plt.title('IQD Equatorial') 
    328     plt.subplot(224) 
    329     plt.plot(q,IQDe/Sqp2-1) 
    330     plt.tight_layout() 
    331     plt.show() 
    332 #comparing monodisperse Beta approximation to version without 
    333     q=np.linspace(1e-5,1) 
    334     F1, F2, IQM, IQSM, IQBM, SQ, SQ_EFF = ellipsoid_Theta(q,20,10,2,1,0.15,12.59921049894873) 
    335     plt.subplot(321) 
    336     plt.loglog(q,IQBM,'g',label='IQBM') 
    337     plt.loglog(q,IQSM,'r',label= 'IQSM') 
    338     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    339     plt.subplot(323) 
    340     plt.plot(q,IQBM,'g',label='IQBM') 
    341     plt.plot(q,IQSM,'r',label= 'IQSM') 
    342     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    343     plt.subplot(325) 
    344     plt.plot(q,IQBM/IQSM,'r',label= 'IQBM/IQSM') 
    345     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    346     plt.tight_layout() 
    347     plt.show() 
    348 #comparing polydispersed Beta approximation to version without 
    349     IQD, IQSD, IQBD, SQ, SQ_EFF = ellipsoid_pe(q, 20,10, 2,1, 0.15,12.59921049894873) 
    350     plt.subplot(321) 
    351     plt.loglog(q, SQ) 
    352     plt.loglog(q, SQ_EFF,label='SQ, SQ_EFF') 
    353     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    354     plt.subplot(323) 
    355     plt.loglog(q,IQBD) 
    356     plt.loglog(q,IQSD,label='IQBD,IQSD') 
    357     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    358     plt.subplot(325) 
    359     plt.plot(q,IQBD) 
    360     plt.plot(q,IQSD,label='IQBD,IQSD') 
    361     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    362     plt.tight_layout() 
    363     plt.show() 
    364  
    365 # YUN 
    366 if (ELLIPSOID == True) & (GAUSSIAN == True) & (YUN == True): 
    367     #Yun's data for Beta Approximation 
    368     data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T    
    369     print('\n\n ELLIPSOID COMPARISON WITH YUN WITHOUT V') 
    370     Q=data[0] 
    371     F1,F2,IQM,IQSM,IQBM,SQ,SQ_EFF = ellipsoid_Theta(Q,20,10,2,1,0.15,12.59921049894873) 
    372     plt.subplot(211) 
    373     plt.loglog(Q,0.15*data[3]*data[5]*ellipsoid_volume(20,10)*1e-4) 
    374     plt.loglog(Q,IQSM,'-g',label='IQSM') 
    375     plt.loglog(data[0], data[3]*ellipsoid_volume(20,10)*1e-4,'-k') 
    376     plt.loglog(Q,IQM,'r',label = 'IQM') 
    377     plt.ylim(1e-5,4) 
    378     plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 
    379     plt.show() 
    380     plt.subplot(221) 
    381     plt.loglog(data[0],SQ) 
    382     plt.loglog(data[0],data[5]) 
    383     plt.title('SQ') 
    384     plt.subplot(222) 
    385     plt.semilogx(data[0],SQ/data[5]-1) 
    386     plt.subplot(223) 
    387     plt.title('SQ_EFF') 
    388     plt.loglog(data[0],SQ_EFF) 
    389     plt.loglog(data[0],data[6]) 
    390     plt.subplot(224) 
    391     plt.semilogx(data[0],(SQ_EFF/data[6])-1,label='Sq effective ratio') 
    392     plt.tight_layout() 
    393     plt.show() 
    394  
    395 #SASFIT 
    396 if ELLIPSOID and GAUSSIAN and SASFIT: 
    397     print('\n\n ELLIPSOID COMPARISON WITH SASFIT WITHOUT V') 
    398     #Rp=20,Re=10,eta_core=4,eta_solv=1 
    399     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQM.txt'),dtype=str,delimiter=';').T 
    400     sasqdata=map(float,sasdata[0]) 
    401     sasvaluedata=map(float,sasdata[1]) 
    402     plt.subplot(321) 
    403     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0)[0],'b') 
    404     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata),'g') 
    405     plt.title('IQM') 
    406     plt.subplot(322) 
    407     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0)[0]-1) 
    408     plt.tight_layout() 
    409     plt.show() 
    410     print('========================================================================================') 
    411     print('========================================================================================') 
    412 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15 
    413     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq.txt'),dtype=str,delimiter=';').T 
    414     sasqdata=map(float,sasdata[0]) 
    415     sasvaluedata=map(float,sasdata[1]) 
    416     plt.subplot(321) 
    417     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]) 
    418     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    419     plt.title('SQ poly 10% 0.15') 
    420     plt.subplot(322) 
    421     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]-1) 
    422 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15 
    423     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq2.txt'),dtype=str,delimiter=';').T 
    424     sasqdata=map(float,sasdata[0]) 
    425     sasvaluedata=map(float,sasdata[1]) 
    426     plt.subplot(323) 
    427     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]) 
    428     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    429     plt.title('SQ poly 30% .15') 
    430     plt.subplot(324) 
    431     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]-1) 
    432 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15 
    433     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq3.txt'),dtype=str,delimiter=';').T 
    434     sasqdata=map(float,sasdata[0]) 
    435     sasvaluedata=map(float,sasdata[1]) 
    436     plt.subplot(325) 
    437     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]) 
    438     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    439     plt.title('SQ poly 60% .15') 
    440     plt.subplot(326) 
    441     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]-1) 
    442     plt.tight_layout() 
    443     plt.show()  
    444     print('========================================================================================') 
    445     print('========================================================================================') 
    446 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3 
    447     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq4.txt'),dtype=str,delimiter=';').T 
    448     sasqdata=map(float,sasdata[0]) 
    449     sasvaluedata=map(float,sasdata[1]) 
    450     plt.subplot(321) 
    451     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]) 
    452     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    453     plt.title('SQ poly 10% .3') 
    454     plt.subplot(322) 
    455     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]-1) 
    456 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3 
    457     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq5.txt'),dtype=str,delimiter=';').T 
    458     sasqdata=map(float,sasdata[0]) 
    459     sasvaluedata=map(float,sasdata[1]) 
    460     plt.subplot(323) 
    461     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]) 
    462     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    463     plt.title('SQ poly 30% .3') 
    464     plt.subplot(324) 
    465     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]-1) 
    466 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3 
    467     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq6.txt'),dtype=str,delimiter=';').T 
    468     sasqdata=map(float,sasdata[0]) 
    469     sasvaluedata=map(float,sasdata[1]) 
    470     plt.subplot(325) 
    471     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]) 
    472     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    473     plt.title('SQ poly 60% .3') 
    474     plt.subplot(326) 
    475     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]-1) 
    476     plt.tight_layout() 
    477     plt.show() 
    478     print('========================================================================================') 
    479     print('========================================================================================') 
    480 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6 
    481     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq7.txt'),dtype=str,delimiter=';').T 
    482     sasqdata=map(float,sasdata[0]) 
    483     sasvaluedata=map(float,sasdata[1]) 
    484     plt.subplot(321) 
    485     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]) 
    486     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    487     plt.title('SQ poly 10% .6') 
    488     plt.subplot(322) 
    489     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]-1) 
    490 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6 
    491     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq8.txt'),dtype=str,delimiter=';').T 
    492     sasqdata=map(float,sasdata[0]) 
    493     sasvaluedata=map(float,sasdata[1]) 
    494     plt.subplot(323) 
    495     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]) 
    496     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    497     plt.title('SQ poly 30% .6') 
    498     plt.subplot(324) 
    499     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]-1) 
    500 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6 
    501     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq9.txt'),dtype=str,delimiter=';').T 
    502     sasqdata=map(float,sasdata[0]) 
    503     sasvaluedata=map(float,sasdata[1]) 
    504     plt.subplot(325) 
    505     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]) 
    506     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    507     plt.title('SQ poly 60% .6') 
    508     plt.subplot(326) 
    509     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]-1) 
    510     plt.tight_layout() 
    511     plt.show() 
    512     print('========================================================================================') 
    513     print('========================================================================================') 
    514 #checks beta structre factor 
    515 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15 
    516     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt'),dtype=str,delimiter=';').T 
    517     sasqdata=map(float,sasdata[0]) 
    518     sasvaluedata=map(float,sasdata[1]) 
    519     plt.subplot(321) 
    520     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]) 
    521     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    522     plt.title('SQ_EFF poly 10% 0.15') 
    523     plt.subplot(322) 
    524     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]-1) 
    525 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15 
    526     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt'),dtype=str,delimiter=';').T 
    527     sasqdata=map(float,sasdata[0]) 
    528     sasvaluedata=map(float,sasdata[1]) 
    529     plt.subplot(323) 
    530     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]) 
    531     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    532     plt.title('SQ_EFF poly 30% .15') 
    533     plt.subplot(324) 
    534     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]-1) 
    535 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15 
    536     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt'),dtype=str,delimiter=';').T 
    537     sasqdata=map(float,sasdata[0]) 
    538     sasvaluedata=map(float,sasdata[1]) 
    539     plt.subplot(325) 
    540     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]) 
    541     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    542     plt.title('SQ_EFF poly 60% .15') 
    543     plt.subplot(326) 
    544     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]-1) 
    545     plt.tight_layout() 
    546     plt.show()  
    547     print('========================================================================================') 
    548     print('========================================================================================') 
    549 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3 
    550     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt'),dtype=str,delimiter=';').T 
    551     sasqdata=map(float,sasdata[0]) 
    552     sasvaluedata=map(float,sasdata[1]) 
    553     plt.subplot(321) 
    554     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]) 
    555     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    556     plt.title('SQ_EFF poly 10% .3') 
    557     plt.subplot(322) 
    558     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]-1) 
    559 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3 
    560     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt'),dtype=str,delimiter=';').T 
    561     sasqdata=map(float,sasdata[0]) 
    562     sasvaluedata=map(float,sasdata[1]) 
    563     plt.subplot(323) 
    564     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]) 
    565     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    566     plt.title('SQ_EFF poly 30% .3') 
    567     plt.subplot(324) 
    568     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]-1) 
    569 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3 
    570     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt'),dtype=str,delimiter=';').T 
    571     sasqdata=map(float,sasdata[0]) 
    572     sasvaluedata=map(float,sasdata[1]) 
    573     plt.subplot(325) 
    574     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]) 
    575     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    576     plt.title('SQ_EFF poly 60% .3') 
    577     plt.subplot(326) 
    578     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]-1) 
    579     plt.tight_layout() 
    580     plt.show() 
    581     print('========================================================================================') 
    582     print('========================================================================================') 
    583 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6 
    584     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt'),dtype=str,delimiter=';').T 
    585     sasqdata=map(float,sasdata[0]) 
    586     sasvaluedata=map(float,sasdata[1]) 
    587     plt.subplot(321) 
    588     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]) 
    589     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    590     plt.title('SQ_EFF poly 10% .6') 
    591     plt.subplot(322) 
    592     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]-1) 
    593 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6 
    594     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt'),dtype=str,delimiter=';').T 
    595     sasqdata=map(float,sasdata[0]) 
    596     sasvaluedata=map(float,sasdata[1]) 
    597     plt.subplot(323) 
    598     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]) 
    599     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    600     plt.title('SQ_EFF poly 30% .6') 
    601     plt.subplot(324) 
    602     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]-1) 
    603 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6 
    604     sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt'),dtype=str,delimiter=';').T 
    605     sasqdata=map(float,sasdata[0]) 
    606     sasvaluedata=map(float,sasdata[1]) 
    607     plt.subplot(325) 
    608     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]) 
    609     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    610     plt.title('SQ_EFF poly 60% .6') 
    611     plt.subplot(326) 
    612     plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]-1) 
    613     plt.tight_layout() 
    614     plt.show() 
    615     print('========================================================================================') 
    616     print('========================================================================================') 
    617 #checks polydispersion for single parameter and varying pd with sasfit 
    618 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    619     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD.txt'),dtype=str,delimiter=';').T 
    620     sasqdata=map(float,sasdata[0]) 
    621     sasvaluedata=map(float,sasdata[1]) 
    622     plt.subplot(221) 
    623     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.1,radius_equatorial_pd=0)[0]) 
    624     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    625     plt.title('IQD polar 10%') 
    626     plt.subplot(222) 
    627     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.1,radius_equatorial_pd=0)[0]-1) 
    628 #N=1,s=1,X0=10,distr equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    629     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD2.txt'),dtype=str,delimiter=';').T 
    630     sasqdata=map(float,sasdata[0]) 
    631     sasvaluedata=map(float,sasdata[1]) 
    632     plt.subplot(223) 
    633     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.1)[0]) 
    634     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    635     plt.title('IQD equatorial 10%') 
    636     plt.subplot(224) 
    637     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.1)[0]-1) 
    638     plt.tight_layout() 
    639     plt.show() 
    640     print('========================================================================================') 
    641     print('========================================================================================') 
    642     #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    643     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD3.txt'),dtype=str,delimiter=';').T 
    644     sasqdata=map(float,sasdata[0]) 
    645     sasvaluedata=map(float,sasdata[1]) 
    646     plt.subplot(221) 
    647     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.3,radius_equatorial_pd=0)[0]) 
    648     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    649     plt.title('IQD polar 30%') 
    650     plt.subplot(222) 
    651     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.3,radius_equatorial_pd=0)[0]-1) 
    652 #N=1,s=3,X0=10,distr equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    653     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD4.txt'),dtype=str,delimiter=';').T 
    654     sasqdata=map(float,sasdata[0]) 
    655     sasvaluedata=map(float,sasdata[1]) 
    656     plt.subplot(223) 
    657     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.3)[0]) 
    658     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    659     plt.title('IQD equatorial 30%') 
    660     plt.subplot(224) 
    661     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.3)[0]-1) 
    662     plt.tight_layout() 
    663     plt.show() 
    664     print('========================================================================================') 
    665     print('========================================================================================') 
    666     #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    667     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD5.txt'),dtype=str,delimiter=';').T 
    668     sasqdata=map(float,sasdata[0]) 
    669     sasvaluedata=map(float,sasdata[1]) 
    670     plt.subplot(221) 
    671     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.6,radius_equatorial_pd=0)[0]) 
    672     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    673     plt.title('IQD polar 60%') 
    674     plt.subplot(222) 
    675     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.6,radius_equatorial_pd=0)[0]-1) 
    676 #N=1,s=6,X0=10,distr equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
    677     sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD6.txt'),dtype=str,delimiter=';').T 
    678     sasqdata=map(float,sasdata[0]) 
    679     sasvaluedata=map(float,sasdata[1]) 
    680     plt.subplot(223) 
    681     plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.6)[0]) 
    682     plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)) 
    683     plt.title('IQD equatorial 60%') 
    684     plt.subplot(224) 
    685     plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.6)[0]-1) 
    686     plt.tight_layout() 
    687     plt.show() 
    688     print('========================================================================================') 
    689     print('========================================================================================') 
    690  
    691 # TEST FOR RICHARD 
     407    model = 'ellipsoid' 
     408    pars = dict( 
     409        radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1, 
     410        background=0, 
     411        radius_polar_pd=.1, radius_polar_pd_type=pd_type, 
     412        radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type, 
     413        volfraction=0.15, 
     414        radius_effective=270.7543927018, 
     415        #radius_effective=12.59921049894873, 
     416        ) 
     417    target = sasmodels_theory(q, model, beta_mode=1, **pars) 
     418    actual = ellipsoid_pe(q, norm='sasview', **pars) 
     419    print(actual) 
     420    title = " ".join(("sasmodels", model, pd_type)) 
     421    compare(title, target, actual) 
     422COMPARISON[('sasview','ellipsoid','gaussian')] = lambda: compare_sasview_ellipsoid(pd_type='gaussian') 
     423COMPARISON[('sasview','ellipsoid','schulz')] = lambda: compare_sasview_ellipsoid(pd_type='schulz') 
     424 
     425def compare_yun_ellipsoid_mono(): 
     426    pars = { 
     427        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian', 
     428        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian', 
     429        'sld': 2, 'sld_solvent': 1, 
     430        'volfraction': 0.15, 
     431        # Yun uses radius for same volume sphere for effective radius 
     432        # whereas sasview uses the average curvature. 
     433        'radius_effective': 12.59921049894873, 
     434    } 
     435 
     436    data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T 
     437    Q = data[0] 
     438    F1 = data[1] 
     439    P = data[3]*pars['volfraction'] 
     440    S = data[5] 
     441    Seff = data[6] 
     442    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) 
     443    actual = ellipsoid_pe(Q, norm='yun', **pars) 
     444    title = " ".join(("yun", "ellipsoid", "no pd")) 
     445    #compare(title, target, actual, fields="P S I Seff Ibeta") 
     446    compare(title, target, actual) 
     447COMPARISON[('yun','ellipsoid','gaussian')] = compare_yun_ellipsoid_mono 
     448COMPARISON[('yun','ellipsoid','schulz')] = compare_yun_ellipsoid_mono 
     449 
     450def compare_yun_sphere_gauss(): 
     451    # Note: yun uses gauss limits from R0/10 to R0 + 5 sigma steps sigma/100 
     452    # With pd = 0.1, that's 14 sigma and 1400 points. 
     453    pars = { 
     454        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian', 
     455        'sld': 6, 'sld_solvent': 0, 
     456        'volfraction': 0.1, 
     457    } 
     458 
     459    data = np.loadtxt(data_file('testPolydisperseGaussianSphere.dat'),skiprows=2).T 
     460    Q = data[0] 
     461    F1 = data[1] 
     462    P = data[3] 
     463    S = data[5] 
     464    Seff = data[6] 
     465    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) 
     466    actual = sphere_r(Q, norm='yun', **pars) 
     467    title = " ".join(("yun", "sphere", "10% dispersion 10% Vf")) 
     468    compare(title, target, actual) 
     469    data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'),skiprows=2).T 
     470    pars.update(radius_pd=0.15) 
     471    Q = data[0] 
     472    F1 = data[1] 
     473    P = data[3] 
     474    S = data[5] 
     475    Seff = data[6] 
     476    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) 
     477    actual = sphere_r(Q, norm='yun', **pars) 
     478    title = " ".join(("yun", "sphere", "15% dispersion 10% Vf")) 
     479    compare(title, target, actual) 
     480COMPARISON[('yun','sphere','gaussian')] = compare_yun_sphere_gauss 
     481 
     482 
     483def compare_sasfit_sphere_gauss(): 
     484    #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3 
     485    pars = { 
     486        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian', 
     487        'sld': 4, 'sld_solvent': 1, 
     488        'volfraction': 0.3, 
     489    } 
     490 
     491    Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt')) 
     492    Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt')) 
     493    Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt')) 
     494    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt')) 
     495    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt')) 
     496    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) 
     497    actual = sphere_r(Q, norm="sasfit", **pars) 
     498    title = " ".join(("sasfit", "sphere", "pd=10% gaussian")) 
     499    compare(title, target, actual) 
     500    #compare(title, target, actual, fields="P") 
     501COMPARISON[('sasfit','sphere','gaussian')] = compare_sasfit_sphere_gauss 
     502 
     503def compare_sasfit_sphere_schulz(): 
    692504    #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1 
    693505    #We have scaled the output from sasfit by 1e-4*volume*volfraction 
    694506    #0.10050378152592121 
    695 if (SPHERE == True) & (SCHULZ == True) & (SASVIEW == True) & (SASFIT == True): 
    696     print('                   *****SPHERE MODEL*****') 
    697     sasdata = np.loadtxt(data_file('richard_test.txt'),dtype=str,delimiter=';').T 
    698     sasqdata=map(float,sasdata[0]) 
    699     model=build_model('sphere',sasqdata) 
    700     Iq = model(radius=20,radius_pd=0.1,radius_pd_n=80,sld=4,sld_solvent=1,background=0,radius_pd_type='schulz',radius_pd_nsigma=8) 
    701     model2=build_model('sphere@hardsphere',sasqdata) 
    702     IQS = model2(radius=20,radius_pd=0.1,radius_pd_n=80,sld=4,sld_solvent=1,background=0,radius_pd_type='schulz',radius_pd_nsigma=8,volfraction=0.3) 
    703      
    704     sasvaluedata=map(float,sasdata[1]) 
    705     sasdata2 = np.loadtxt(data_file('richard_test2.txt'),dtype=str,delimiter=';').T 
    706     sasqdata2=map(float,sasdata2[0]) 
    707     sasvaluedata2=map(float,sasdata2[1]) 
    708     sasdata3 = np.loadtxt(data_file('richard_test3.txt'),dtype=str,delimiter=';').T 
    709     sasqdata3=map(float,sasdata3[0]) 
    710     sasvaluedata3=map(float,sasdata3[1]) 
    711     IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(np.array(sasqdata),20,4,1,.3,0.1,distribution='schulz') 
    712     plt.grid(True) 
    713     plt.title('IQD(blue) vs. SASVIEW(red) vs. SASFIT(green)') 
    714     plt.loglog(sasqdata,Iq,'b') 
    715     plt.loglog(sasqdata,IQD,'r') 
    716     plt.loglog(sasqdata,1e-4/(4./3.*pi*20**3)*np.array(sasvaluedata),'g') 
    717     plt.show() 
    718     plt.grid(True) 
    719     plt.title('IQSD(blue) vs. SASVIEW(red) vs. SASFIT(green)') 
    720     plt.loglog(sasqdata,IQSD,'b') 
    721     plt.loglog(sasqdata,IQS,'r') 
    722     plt.loglog(sasqdata,1e-4/(4./3.*pi*20**3)*0.3*np.array(sasvaluedata2),'g') 
    723     plt.show() 
    724     plt.grid(True) 
    725     plt.title('IQBD(blue) vs. SASFIT(green)') 
    726     plt.loglog(sasqdata,IQBD,'b') 
    727     plt.loglog(sasqdata,1e-4/(4./3.*pi*20**3)*0.3*np.array(sasvaluedata3),'g') 
    728     plt.show() 
    729  
    730 # TEST FOR RICHARD 
    731 if (ELLIPSOID == True) & (SCHULZ == True) & (SASVIEW == True) & (SASFIT == True):    
     507    pars = { 
     508        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz', 
     509        'sld': 4, 'sld_solvent': 1, 
     510        'volfraction': 0.3, 
     511    } 
     512 
     513    Q, IQ = load_sasfit(data_file('richard_test.txt')) 
     514    Q, IQSD = load_sasfit(data_file('richard_test2.txt')) 
     515    Q, IQBD = load_sasfit(data_file('richard_test3.txt')) 
     516    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
     517    actual = sphere_r(Q, norm="sasfit", **pars) 
     518    title = " ".join(("sasfit", "sphere", "pd=10% schulz")) 
     519    compare(title, target, actual) 
     520COMPARISON[('sasfit','sphere','schulz')] = compare_sasfit_sphere_schulz 
     521 
     522def compare_sasfit_ellipsoid_schulz(): 
    732523    #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1 
    733          #Effective radius =13.1353356684 
    734         #We have scaled the output from sasfit by 1e-4*volume*volfraction 
    735         #0.10050378152592121 
    736     print('                   *****ELLIPSOID MODEL*****') 
    737     sasdata = np.loadtxt(data_file('richard_test4.txt'),dtype=str,delimiter=';').T 
    738     sasqdata=map(float,sasdata[0]) 
    739     model=build_model('ellipsoid',sasqdata) 
    740     Iq = model(radius_polar=20,radius_polar_pd=0.1,radius_polar_pd_n=80, radius_equatorial=10, radius_equatorial_pd=0, sld=4,sld_solvent=1,background=0,radius_polar_pd_type='schulz',radius_polar_pd_nsigma=8) 
    741     model2=build_model('ellipsoid@hardsphere',sasqdata) 
    742     IQS =  model2(radius_polar=20,radius_equatorial=10,sld=4,sld_solvent=1,background=0,radius_polar_pd=.1,\ 
    743                      radius_polar_pd_n=80,radius_equatorial_pd=0,radius_polar_pd_type='schulz',radius_polar_pd_nsigma=8,volfraction=0.3) 
    744     sasvaluedata=map(float,sasdata[1]) 
    745     sasdata2 = np.loadtxt(data_file('richard_test5.txt'),dtype=str,delimiter=';').T 
    746     sasqdata2=map(float,sasdata2[0]) 
    747     sasvaluedata2=map(float,sasdata2[1]) 
    748     sasdata3 = np.loadtxt(data_file('richard_test6.txt'),dtype=str,delimiter=';').T 
    749     sasqdata3=map(float,sasdata3[0]) 
    750     sasvaluedata3=map(float,sasdata3[1]) 
    751     IQD, IQSD, IQBD, SQ, SQ_EFF = ellipsoid_pe(np.array(sasqdata),20,10,4,1,.3,radius_polar_pd=0.1,radius_equatorial_pd=0,distribution='schulz')     
    752     plt.grid(True) 
    753     plt.title('IQD(blue) vs. SASVIEW(red) vs. SASFIT(green)') 
    754     plt.loglog(sasqdata,Iq,'b') 
    755     plt.loglog(sasqdata,IQD,'r') 
    756     plt.loglog(sasqdata,1e-4/(4./3.*pi*20*10**2)*np.array(sasvaluedata),'g') 
    757     plt.show() 
    758     plt.grid(True) 
    759     plt.title('IQSD(blue) vs. SASVIEW(red) vs. SASFIT(green)') 
    760     plt.loglog(sasqdata,IQSD,'b') 
    761     plt.loglog(sasqdata,IQS,'r') 
    762     plt.loglog(sasqdata,1e-4/(4./3.*pi*20*10**2)*0.3*np.array(sasvaluedata2),'g') 
    763     plt.show() 
    764     plt.grid(True) 
    765     plt.title('IQBD(blue) vs. SASFIT(green)') 
    766     plt.loglog(sasqdata,IQBD,'b') 
    767     plt.loglog(sasqdata,1e-4/(4./3.*pi*20*10**2)*0.3*np.array(sasvaluedata3),'g') 
    768     plt.show()     
    769      
    770 # SASVIEW GAUSS 
    771 if (SPHERE == True) & (GAUSSIAN == True) & (SASVIEW == True): 
    772     q = np.logspace(-5, 0, 200) 
    773     IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(q,20,4,1,.3) 
    774     model = build_model("sphere", q) 
    775     Sq = model(radius=20,radius_pd=0.1,radius_pd_n=35,sld=4,sld_solvent=1,background=0) 
    776     model2 = build_model("sphere@hardsphere", q) 
    777     S=build_model("hardsphere",q)(radius_effective=20,volfraction=.3,background=0) 
    778     Sq2 = model2(radius=20,radius_pd=0.1,radius_pd_n=35,sld=4,sld_solvent=1,background=0,volfraction=.3) 
    779     print('\n\n SPHERE COMPARISON WITH SASVIEW WITHOUT V') 
    780     plt.subplot(221) 
    781     plt.title('IQD') 
    782     plt.loglog(q, IQD, '-b') 
    783     plt.loglog(q, Sq, '-r') 
    784     plt.subplot(222) 
    785     plt.semilogx(q, Sq/IQD-1, '-g') 
    786     plt.tight_layout() 
    787     plt.show() 
    788     plt.subplot(221) 
    789     plt.title('SQ') 
    790     plt.plot(q, SQ, '-r') 
    791     plt.plot(q,S,'-k') 
    792     plt.subplot(222) 
    793     plt.plot(SQ/S-1) 
    794     plt.tight_layout() 
    795     plt.show() 
    796     plt.subplot(221) 
    797     plt.title('IQSD') 
    798     plt.loglog(q, IQSD, '-b') 
    799     plt.loglog(q, Sq2, '-r',label='IQSD') 
    800     plt.subplot(222) 
    801     plt.semilogx(q, Sq2/IQSD-1, '-g',label='IQSD ratio') 
    802     plt.tight_layout() 
    803     plt.show() 
    804     IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(q,20,4,1,0.15) 
    805     plt.subplot(211) 
    806     plt.title('SQ vs SQ_EFF') 
    807     plt.plot(q, SQ) 
    808     plt.plot(q, SQ_EFF) 
    809     plt.tight_layout() 
    810     plt.show()  
    811     plt.subplot(221) 
    812     plt.title('IQSD vs IQBD') 
    813     plt.loglog(q,IQBD) 
    814     plt.loglog(q,IQSD,label='IQBD,IQSD') 
    815     plt.subplot(222) 
    816     plt.plot(q,IQBD) 
    817     plt.plot(q,IQSD,) 
    818     plt.tight_layout() 
    819     plt.show()  
    820  
    821 # SASFIT GAUSS 
    822 if (SPHERE == True) & (GAUSSIAN == True) & (SASFIT == True): 
    823     print('\n\n SPHERE COMPARISON WITH SASFIT WITHOUT V') 
    824 #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3 
    825     sasdata = np.loadtxt(data_file('sasfit_sphere_IQD.txt'),dtype=str,delimiter=';').T 
    826     sasqdata=map(float,sasdata[0]) 
    827     sasvaluedata=map(float,sasdata[1]) 
    828     IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(np.array(sasqdata),20,4,1,.3) 
    829     plt.subplot(221) 
    830     plt.title('IQD') 
    831     plt.loglog(sasqdata,IQD ) 
    832     plt.loglog(sasqdata,1e-4/(4./3*pi*20**3)*np.array(sasvaluedata)) 
    833     plt.subplot(222) 
    834     plt.semilogx(sasqdata,1e-4/(4./3*pi*20**3)*np.array(sasvaluedata)/IQD-1) 
    835     plt.tight_layout() 
    836     plt.show() 
    837     sasdata = np.loadtxt(data_file('sasfit_sphere_IQSD.txt'),dtype=str,delimiter=';').T 
    838     sasqdata=map(float,sasdata[0]) 
    839     sasvaluedata=map(float,sasdata[1]) 
    840     plt.subplot(221) 
    841     plt.title('IQSD') 
    842     plt.loglog(sasqdata, IQSD) 
    843     plt.loglog(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata)) 
    844     plt.subplot(222) 
    845     plt.semilogx(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata)/IQSD-1) 
    846     plt.tight_layout() 
    847     plt.show() 
    848     sasdata = np.loadtxt(data_file('sasfit_sphere_IQBD.txt'),dtype=str,delimiter=';').T 
    849     sasqdata=map(float,sasdata[0]) 
    850     sasvaluedata=map(float,sasdata[1]) 
    851     plt.subplot(221) 
    852     plt.title('IQBD') 
    853     plt.loglog(sasqdata,IQBD) 
    854     plt.loglog(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata)) 
    855     plt.subplot(222) 
    856     plt.semilogx(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata)/IQBD-1) 
    857     plt.tight_layout() 
    858     plt.show() 
    859     sasdata = np.loadtxt(data_file('sasfit_polydisperse_sphere_sq.txt'),dtype=str,delimiter=';').T 
    860     sasqdata=map(float,sasdata[0]) 
    861     sasvaluedata=map(float,sasdata[1]) 
    862     plt.subplot(221) 
    863     plt.title('SQ') 
    864     plt.loglog(sasqdata, SQ) 
    865     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    866     plt.subplot(222) 
    867     plt.semilogx(sasqdata,np.array(sasvaluedata)/SQ-1) 
    868     plt.tight_layout() 
    869     plt.show() 
    870     sasdata = np.loadtxt(data_file('sasfit_polydisperse_sphere_sqeff.txt'),dtype=str,delimiter=';').T 
    871     sasqdata=map(float,sasdata[0]) 
    872     sasvaluedata=map(float,sasdata[1]) 
    873     plt.subplot(221) 
    874     plt.title('SQ_EFF') 
    875     plt.loglog(sasqdata,SQ_EFF) 
    876     plt.loglog(sasqdata,np.array(sasvaluedata)) 
    877     plt.subplot(222) 
    878     plt.semilogx(sasqdata,np.array(sasvaluedata)/SQ_EFF-1) 
    879     plt.tight_layout() 
    880     plt.show() 
     524    #Effective radius =13.1353356684 
     525    #We have scaled the output from sasfit by 1e-4*volume*volfraction 
     526    #0.10050378152592121 
     527    pars = { 
     528        'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz', 
     529        'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz', 
     530        'sld': 4, 'sld_solvent': 1, 
     531        'volfraction': 0.3, 'radius_effective': 13.1353356684, 
     532    } 
     533 
     534    Q, IQ = load_sasfit(data_file('richard_test4.txt')) 
     535    Q, IQSD = load_sasfit(data_file('richard_test5.txt')) 
     536    Q, IQBD = load_sasfit(data_file('richard_test6.txt')) 
     537    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
     538    actual = ellipsoid_pe(Q, norm="sasfit", **pars) 
     539    title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz")) 
     540    compare(title, target, actual) 
     541COMPARISON[('sasfit','ellipsoid','schulz')] = compare_sasfit_ellipsoid_schulz 
     542 
     543 
     544def compare_sasfit_ellipsoid_gaussian(): 
     545    pars = { 
     546        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian', 
     547        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian', 
     548        'sld': 4, 'sld_solvent': 1, 
     549        'volfraction': 0, 'radius_effective': None, 
     550    } 
     551 
     552    #Rp=20,Re=10,eta_core=4,eta_solv=1 
     553    Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt')) 
     554    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None) 
     555    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     556    target = Theory(Q=Q, P=PQ0) 
     557    compare("sasfit ellipsoid no poly", target, actual); plt.show() 
     558 
     559    #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     560    Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt')) 
     561    pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None) 
     562    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     563    target = Theory(Q=Q, P=PQ_Rp10) 
     564    compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show() 
     565    #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     566    Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt')) 
     567    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None) 
     568    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     569    target = Theory(Q=Q, P=PQ_Re10) 
     570    compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show() 
     571    #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     572    Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt')) 
     573    pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None) 
     574    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     575    target = Theory(Q=Q, P=PQ_Rp30) 
     576    compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show() 
     577    #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     578    Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt')) 
     579    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None) 
     580    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     581    target = Theory(Q=Q, P=PQ_Re30) 
     582    compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show() 
     583    #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     584    Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt')) 
     585    pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None) 
     586    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     587    target = Theory(Q=Q, P=PQ_Rp60) 
     588    compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show() 
     589    #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 
     590    Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt')) 
     591    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None) 
     592    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     593    target = Theory(Q=Q, P=PQ_Re60) 
     594    compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show() 
     595 
     596    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15 
     597    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt')) 
     598    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt')) 
     599    pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) 
     600    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     601    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     602    compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show() 
     603    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15 
     604    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt')) 
     605    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt')) 
     606    pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) 
     607    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     608    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     609    compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show() 
     610    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15 
     611    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt')) 
     612    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt')) 
     613    pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) 
     614    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     615    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     616    compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show() 
     617 
     618    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3 
     619    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt')) 
     620    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt')) 
     621    pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) 
     622    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     623    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     624    compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show() 
     625    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3 
     626    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt')) 
     627    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt')) 
     628    pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) 
     629    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     630    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     631    compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show() 
     632    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3 
     633    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt')) 
     634    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt')) 
     635    pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) 
     636    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     637    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     638    compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show() 
     639 
     640    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6 
     641    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt')) 
     642    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt')) 
     643    pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) 
     644    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     645    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     646    compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show() 
     647    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6 
     648    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt')) 
     649    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt')) 
     650    pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) 
     651    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     652    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     653    compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show() 
     654    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6 
     655    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt')) 
     656    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt')) 
     657    pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) 
     658    actual = ellipsoid_pe(Q, norm='sasfit', **pars) 
     659    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 
     660    compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show() 
     661COMPARISON[('sasfit','ellipsoid','gaussian')] = compare_sasfit_ellipsoid_gaussian 
     662 
     663def main(): 
     664    key = tuple(sys.argv[1:]) 
     665    if key not in COMPARISON: 
     666        print("usage: sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid] [gaussian|schulz]") 
     667        return 
     668    comparison = COMPARISON[key] 
     669    comparison() 
     670 
     671if __name__ == "__main__": 
     672    main() 
  • sasmodels/core.py

    r3221de0 r01c8d9e  
    140140    used with functions in generate to build the docs or extract model info. 
    141141    """ 
     142 
    142143    if "+" in model_string: 
    143144        parts = [load_model_info(part) 
     
    205206 
    206207    numpy_dtype, fast, platform = parse_dtype(model_info, dtype, platform) 
    207  
    208208    source = generate.make_source(model_info) 
    209209    if platform == "dll": 
     
    265265    # Assign default platform, overriding ocl with dll if OpenCL is unavailable 
    266266    # If opencl=False OpenCL is switched off 
    267  
    268267    if platform is None: 
    269268        platform = "ocl" 
     
    291290    else: 
    292291        numpy_dtype = np.dtype(dtype) 
    293  
    294292    # Make sure that the type is supported by opencl, otherwise use dll 
    295293    if platform == "ocl": 
  • sasmodels/data.py

    r65fbf7c r01c8d9e  
    329329    else: 
    330330        dq = resolution * q 
    331  
    332331    data = Data1D(q, Iq, dx=dq, dy=dIq) 
    333332    data.filename = "fake data" 
     
    523522                  else 'log') 
    524523        plt.xscale(xscale) 
     524         
    525525        plt.xlabel("$q$/A$^{-1}$") 
    526526        plt.yscale(yscale) 
  • sasmodels/generate.py

    rd86f0fc r01c8d9e  
    672672 
    673673# type in IQXY pattern could be single, float, double, long double, ... 
    674 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
     674_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]", 
    675675                           flags=re.MULTILINE) 
    676676def find_xy_mode(source): 
     
    701701    return 'qa' 
    702702 
     703# type in IQXY pattern could be single, float, double, long double, ... 
     704_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE) 
     705def has_Fq(source): 
     706    for code in source: 
     707        m = _FQ_PATTERN.search(code) 
     708        if m is not None: 
     709            return True 
     710    return False 
    703711 
    704712def _add_source(source, code, path, lineno=1): 
     
    730738    # dispersion.  Need to be careful that necessary parameters are available 
    731739    # for computing volume even if we allow non-disperse volume parameters. 
    732  
    733740    partable = model_info.parameters 
    734  
    735741    # Load templates and user code 
    736742    kernel_header = load_template('kernel_header.c') 
    737743    kernel_code = load_template('kernel_iq.c') 
    738744    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
    739  
    740745    # Build initial sources 
    741746    source = [] 
     
    743748    for path, code in user_code: 
    744749        _add_source(source, code, path) 
    745  
    746750    if model_info.c_code: 
    747751        _add_source(source, model_info.c_code, model_info.filename, 
     
    789793    source.append("\\\n".join(p.as_definition() 
    790794                              for p in partable.kernel_parameters)) 
    791  
    792795    # Define the function calls 
    793796    if partable.form_volume_parameters: 
     
    800803        call_volume = "#define CALL_VOLUME(v) 1.0" 
    801804    source.append(call_volume) 
    802  
    803805    model_refs = _call_pars("_v.", partable.iq_parameters) 
    804     pars = ",".join(["_q"] + model_refs) 
    805     call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     806    #create varaible BETA to turn on and off beta approximation 
     807    BETA = has_Fq(source) 
     808    if not BETA: 
     809        pars = ",".join(["_q"] + model_refs) 
     810        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     811    else: 
     812        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs) 
     813        call_iq = "#define CALL_IQ(_q, _F1, _F2, _v) Fq(%s)" % pars 
    806814    if xy_mode == 'qabc': 
    807815        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
     
    831839    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    832840               if p.type == 'sld'] 
    833  
    834841    # Fill in definitions for numbers of parameters 
     842    source.append("#define BETA %d" %(1 if BETA else 0)) 
    835843    source.append("#define MAX_PD %s"%partable.max_pd) 
    836844    source.append("#define NUM_PARS %d"%partable.npars) 
     
    839847    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    840848    source.append("#define PROJECTION %d"%PROJECTION) 
    841  
    842849    # TODO: allow mixed python/opencl kernels? 
    843  
    844850    ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    845851    dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
     852 
    846853    result = { 
    847854        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    848855        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    849856    } 
    850  
     857    #print(result['dll']) 
    851858    return result 
    852859 
     
    10681075        model_info = make_model_info(kernel_module) 
    10691076        source = make_source(model_info) 
    1070         print(source['dll']) 
     1077        #print(source['dll']) 
    10711078 
    10721079 
  • sasmodels/kernel_iq.c

    r7c35fda r01c8d9e  
    3636//  PROJECTION : equirectangular=1, sinusoidal=2 
    3737//      see explore/jitter.py for definitions. 
     38 
    3839 
    3940#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
     
    270271#endif // _QABC_SECTION 
    271272 
    272  
    273273// ==================== KERNEL CODE ======================== 
    274  
    275274kernel 
    276275void KERNEL_NAME( 
     
    339338  // The code differs slightly between opencl and dll since opencl is only 
    340339  // seeing one q value (stored in the variable "this_result") while the dll 
    341   // version must loop over all q. 
     340  // version must loop over all q 
     341  #if BETA 
     342  double *beta_result = &(result[nq+1]); // = result + nq + 1 
     343  double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     344  #endif 
    342345  #ifdef USE_OPENCL 
    343346    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    344347    double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
     348    #if BETA 
     349      double this_beta_result = (pd_start == 0 ? 0.0 : result[nq + q_index]; 
    345350  #else // !USE_OPENCL 
    346     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     351    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);  
    347352    if (pd_start == 0) { 
    348353      #ifdef USE_OPENMP 
     
    350355      #endif 
    351356      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     357      #if BETA 
     358      for (int q_index=0; q_index < nq; q_index++) beta_result[q_index] = 0.0; 
     359      #endif 
    352360    } 
    353361    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    354362#endif // !USE_OPENCL 
     363 
     364 
     365 
    355366 
    356367 
     
    442453// inner loop and defines the macros that use them. 
    443454 
     455 
    444456#if defined(CALL_IQ) 
    445457  // unoriented 1D 
    446458  double qk; 
    447   #define FETCH_Q() do { qk = q[q_index]; } while (0) 
    448   #define BUILD_ROTATION() do {} while(0) 
    449   #define APPLY_ROTATION() do {} while(0) 
    450   #define CALL_KERNEL() CALL_IQ(qk, local_values.table) 
     459  #if BETA == 0 
     460    #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     461    #define BUILD_ROTATION() do {} while(0) 
     462    #define APPLY_ROTATION() do {} while(0) 
     463    #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 
     464 
     465  // unoriented 1D Beta 
     466  #elif BETA == 1 
     467    double F1, F2; 
     468    #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     469    #define BUILD_ROTATION() do {} while(0) 
     470    #define APPLY_ROTATION() do {} while(0) 
     471    #define CALL_KERNEL() CALL_IQ(qk,F1,F2,local_values.table) 
     472  #endif 
    451473 
    452474#elif defined(CALL_IQ_A) 
     
    647669      pd_norm += weight * CALL_VOLUME(local_values.table); 
    648670      BUILD_ROTATION(); 
    649  
     671#if BETA 
     672    if (weight > cutoff) { 
     673      weight_norm += weight;} 
     674#endif 
    650675#ifndef USE_OPENCL 
    651676      // DLL needs to explicitly loop over the q values. 
     
    693718          } 
    694719        #else  // !MAGNETIC 
    695           const double scattering = CALL_KERNEL(); 
     720          #if defined(CALL_IQ) && BETA == 1 
     721            CALL_KERNEL(); 
     722            const double scatteringF1 = F1; 
     723            const double scatteringF2 = F2; 
     724          #else 
     725            const double scattering = CALL_KERNEL(); 
     726          #endif 
    696727        #endif // !MAGNETIC 
    697728//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 
    698729 
    699730        #ifdef USE_OPENCL 
    700           this_result += weight * scattering; 
     731          #if defined(CALL_IQ)&& BETA == 1 
     732             this_result += weight * scatteringF2; 
     733             this_beta_result += weight * scatteringF1; 
     734            #else 
     735              this_result += weight * scattering; 
     736          #endif 
    701737        #else // !USE_OPENCL 
    702           result[q_index] += weight * scattering; 
     738          #if defined(CALL_IQ)&& BETA == 1 
     739            result[q_index] += weight * scatteringF2; 
     740            beta_result[q_index] += weight * scatteringF1; 
     741            #endif 
     742            #else 
     743            result[q_index] += weight * scattering; 
     744          #endif 
    703745        #endif // !USE_OPENCL 
    704746      } 
    705747    } 
    706748  } 
    707  
    708749// close nested loops 
    709750++step; 
     
    728769  result[q_index] = this_result; 
    729770  if (q_index == 0) result[nq] = pd_norm; 
     771  #if BETA 
     772  beta_result[q_index] = this_beta_result; 
     773  #endif 
     774  if (q_index == 0) result[nq] = pd_norm; 
     775 
    730776//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
    731777#else // !USE_OPENCL 
    732778  result[nq] = pd_norm; 
     779  #if BETA 
     780  result[2*nq+1] = weight_norm; 
     781  #endif 
    733782//printf("res: %g/%g\n", result[0], pd_norm); 
    734783#endif // !USE_OPENCL 
  • sasmodels/kernelcl.py

    rd86f0fc r01c8d9e  
    420420        if self.program is None: 
    421421            compile_program = environment().compile_program 
     422            with open('model.c','w') as fid: 
     423                print(self.source['opencl'], file=fid) 
    422424            timestamp = generate.ocl_timestamp(self.info) 
    423425            self.program = compile_program( 
     
    539541        self.dim = '2d' if q_input.is_2d else '1d' 
    540542        # plus three for the normalization values 
    541         self.result = np.empty(q_input.nq+1, dtype) 
     543        self.result = np.empty(2*q_input.nq+2,dtype) 
    542544 
    543545        # Inputs and outputs for each kernel call 
     
    555557                     else np.float16 if dtype == generate.F16 
    556558                     else np.float32)  # will never get here, so use np.float32 
    557  
    558     def __call__(self, call_details, values, cutoff, magnetic): 
     559    __call__= Iq 
     560 
     561    def Iq(self, call_details, values, cutoff, magnetic): 
    559562        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    560563        context = self.queue.context 
     
    572575        ] 
    573576        #print("Calling OpenCL") 
    574         #call_details.show(values) 
    575         # Call kernel and retrieve results 
     577        call_details.show(values) 
     578        #Call kernel and retrieve results 
    576579        wait_for = None 
    577580        last_nap = time.clock() 
     
    604607        return scale*self.result[:self.q_input.nq] + background 
    605608        # return self.result[:self.q_input.nq] 
     609     #NEEDS TO BE FINISHED FOR OPENCL 
     610     def beta(): 
     611         return 0 
    606612 
    607613    def release(self): 
  • sasmodels/kerneldll.py

    r33969b6 r01c8d9e  
    258258        # Note: if there is a syntax error then compile raises an error 
    259259        # and the source file will not be deleted. 
    260         os.unlink(filename) 
    261         #print("saving compiled file in %r"%filename) 
     260        #os.unlink(filename) 
     261        print("saving compiled file in %r"%filename) 
    262262    return dll 
    263263 
     
    371371    def __init__(self, kernel, model_info, q_input): 
    372372        # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 
     373        #,model_info,q_input) 
    373374        self.kernel = kernel 
    374375        self.info = model_info 
     
    376377        self.dtype = q_input.dtype 
    377378        self.dim = '2d' if q_input.is_2d else '1d' 
    378         self.result = np.empty(q_input.nq+1, q_input.dtype) 
     379        self.result = np.empty(2*q_input.nq+2, q_input.dtype) 
    379380        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    380381                     else np.float64 if self.q_input.dtype == generate.F64 
    381382                     else np.float128) 
    382383 
    383     def __call__(self, call_details, values, cutoff, magnetic): 
     384    def Iq(self, call_details, values, cutoff, magnetic): 
    384385        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    385  
     386        self._call_kernel(call_details, values, cutoff, magnetic) 
     387        #print("returned",self.q_input.q, self.result) 
     388        pd_norm = self.result[self.q_input.nq] 
     389        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
     390        background = values[1] 
     391        #print("scale",scale,background) 
     392        return scale*self.result[:self.q_input.nq] + background 
     393    __call__ = Iq 
     394 
     395    def beta(self, call_details, values, cutoff, magnetic): 
     396        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
     397        self._call_kernel(call_details, values, cutoff, magnetic) 
     398        w_norm = self.result[2*self.q_input.nq + 1] 
     399        pd_norm = self.result[self.q_input.nq] 
     400        if w_norm == 0.: 
     401            w_norm = 1. 
     402        F2 = self.result[:self.q_input.nq]/w_norm 
     403        F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm 
     404        volume_avg = pd_norm/w_norm 
     405        return F1, F2, volume_avg 
     406 
     407    def _call_kernel(self, call_details, values, cutoff, magnetic): 
    386408        kernel = self.kernel[1 if magnetic else 0] 
    387409        args = [ 
     
    390412            None, # pd_stop pd_stride[MAX_PD] 
    391413            call_details.buffer.ctypes.data, # problem 
    392             values.ctypes.data,  #pars 
    393             self.q_input.q.ctypes.data, #q 
     414            values.ctypes.data,  # pars 
     415            self.q_input.q.ctypes.data, # q 
    394416            self.result.ctypes.data,   # results 
    395417            self.real(cutoff), # cutoff 
    396418        ] 
    397         #print("Calling DLL") 
     419        #print(self.beta_result.ctypes.data) 
     420#        print("Calling DLL line 397") 
     421#        print("values", values) 
    398422        #call_details.show(values) 
    399423        step = 100 
     
    403427            kernel(*args) # type: ignore 
    404428 
    405         #print("returned",self.q_input.q, self.result) 
    406         pd_norm = self.result[self.q_input.nq] 
    407         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    408         background = values[1] 
    409         #print("scale",scale,background) 
    410         return scale*self.result[:self.q_input.nq] + background 
    411  
    412429    def release(self): 
    413430        # type: () -> None 
  • sasmodels/modelinfo.py

    r95498a3 r01c8d9e  
    163163    parameter.length = length 
    164164    parameter.length_control = control 
    165  
    166165    return parameter 
    167166 
     
    418417    # scale and background are implicit parameters 
    419418    COMMON = [Parameter(*p) for p in COMMON_PARAMETERS] 
    420  
    421419    def __init__(self, parameters): 
    422420        # type: (List[Parameter]) -> None 
    423421        self.kernel_parameters = parameters 
    424422        self._set_vector_lengths() 
    425  
    426423        self.npars = sum(p.length for p in self.kernel_parameters) 
    427424        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
     
    430427        if self.nmagnetic: 
    431428            self.nvalues += 3 + 3*self.nmagnetic 
    432  
    433429        self.call_parameters = self._get_call_parameters() 
    434430        self.defaults = self._get_defaults() 
     
    444440        self.form_volume_parameters = [p for p in self.kernel_parameters 
    445441                                       if p.type == 'volume'] 
    446  
    447442        # Theta offset 
    448443        offset = 0 
     
    466461        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    467462                                if p.id.startswith('M0:')] 
    468  
    469463        self.pd_1d = set(p.name for p in self.call_parameters 
    470464                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
     
    770764        # Custom sum/multi models 
    771765        return kernel_module.model_info 
     766 
    772767    info = ModelInfo() 
    773768    #print("make parameter table", kernel_module.parameters) 
     
    822817    info.lineno = {} 
    823818    _find_source_lines(info, kernel_module) 
    824  
    825819    return info 
    826820 
  • sasmodels/models/ellipsoid.c

    r108e70e r01c8d9e  
    2020    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
    2121    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
    22  
     22   
    2323    // translate a point in [-1,1] to a point in [0, 1] 
    2424    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     
    3636    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    3737    return 1.0e-4 * s * s * form; 
     38}  
     39 
     40static void 
     41Fq(double q, 
     42    double *F1, 
     43    double *F2, 
     44    double sld, 
     45    double sld_solvent, 
     46    double radius_polar, 
     47    double radius_equatorial) 
     48{ 
     49    // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955) 
     50    //     i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT 
     51    //          = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT 
     52    //          = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT 
     53    // u-substitution of 
     54    //     u = sin, du = cos dT 
     55    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
     56    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
     57    // translate a point in [-1,1] to a point in [0, 1] 
     58    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     59    const double zm = 0.5; 
     60    const double zb = 0.5; 
     61    double total_F2 = 0.0; 
     62    double total_F1 = 0.0; 
     63    for (int i=0;i<GAUSS_N;i++) { 
     64        const double u = GAUSS_Z[i]*zm + zb; 
     65        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
     66        const double f = sas_3j1x_x(q*r); 
     67        total_F2 += GAUSS_W[i] * f * f; 
     68        total_F1 += GAUSS_W[i] * f; 
     69    } 
     70    // translate dx in [-1,1] to dx in [lower,upper] 
     71    const double form_squared_avg = total_F2*zm; 
     72    const double form_avg = total_F1*zm; 
     73    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
     74    *F2 = 1e-4 * s * s * form_squared_avg; 
     75    *F1 = 1e-2 * s * form_avg; 
    3876} 
     77 
     78 
     79 
     80 
     81 
    3982 
    4083static double 
  • sasmodels/models/lib/sphere_form.c

    r925ad6e r01c8d9e  
    1313    return 1.0e-4*square(contrast * fq); 
    1414} 
     15 
  • sasmodels/models/sphere.py

    ref07e95 r01c8d9e  
    6767             ] 
    6868 
    69 source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c"] 
     69source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c", "sphere.c"] 
    7070 
    71 # No volume normalization despite having a volume parameter 
    72 # This should perhaps be volume normalized? 
    73 form_volume = """ 
    74     return sphere_volume(radius); 
    75     """ 
    76  
    77 Iq = """ 
    78     return sphere_form(q, radius, sld, sld_solvent); 
    79     """ 
    8071 
    8172def ER(radius): 
  • sasmodels/product.py

    r2d81cfe r01c8d9e  
    1616import numpy as np  # type: ignore 
    1717 
    18 from .modelinfo import ParameterTable, ModelInfo 
     18from .modelinfo import ParameterTable, ModelInfo, Parameter 
    1919from .kernel import KernelModel, Kernel 
    2020from .details import make_details, dispersion_mesh 
     
    7474    translate_name = dict((old.id, new.id) for old, new 
    7575                          in zip(s_pars.kernel_parameters[1:], s_list)) 
    76     combined_pars = p_pars.kernel_parameters + s_list 
     76    beta_parameter = Parameter("beta_mode", "", 0, [["P*S"],["P*(1+beta*(S-1))"], "", "Structure factor dispersion calculation mode"]) 
     77    combined_pars = p_pars.kernel_parameters + s_list + [beta_parameter] 
    7778    parameters = ParameterTable(combined_pars) 
    7879    parameters.max_pd = p_pars.max_pd + s_pars.max_pd 
     
    151152        #: Structure factor modelling interaction between particles. 
    152153        self.S = S 
     154         
    153155        #: Model precision. This is not really relevant, since it is the 
    154156        #: individual P and S models that control the effective dtype, 
     
    168170        # in opencl; or both in opencl, but one in single precision and the 
    169171        # other in double precision). 
     172         
    170173        p_kernel = self.P.make_kernel(q_vectors) 
    171174        s_kernel = self.S.make_kernel(q_vectors) 
     
    193196        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray 
    194197        p_info, s_info = self.info.composition[1] 
    195  
    196198        # if there are magnetic parameters, they will only be on the 
    197199        # form factor P, not the structure factor S. 
     
    205207        nweights = call_details.num_weights 
    206208        weights = values[nvalues:nvalues + 2*nweights] 
    207  
    208209        # Construct the calling parameters for P. 
    209210        p_npars = p_info.parameters.npars 
     
    218219        p_values.append([0.]*spacer) 
    219220        p_values = np.hstack(p_values).astype(self.p_kernel.dtype) 
    220  
    221221        # Call ER and VR for P since these are needed for S. 
    222222        p_er, p_vr = calc_er_vr(p_info, p_details, p_values) 
    223223        s_vr = (volfrac/p_vr if p_vr != 0. else volfrac) 
    224224        #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr)) 
    225  
    226225        # Construct the calling parameters for S. 
    227226        # The  effective radius is not in the combined parameter list, so 
     
    253252        s_values.append([0.]*spacer) 
    254253        s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 
    255  
     254        # beta mode is the first parameter after the structure factor pars 
     255        beta_index = 2+p_npars+s_npars 
     256        beta_mode = values[beta_index] 
    256257        # Call the kernels 
    257         p_result = self.p_kernel(p_details, p_values, cutoff, magnetic) 
    258         s_result = self.s_kernel(s_details, s_values, cutoff, False) 
    259  
     258        if beta_mode: # beta: 
     259            F1, F2, volume_avg = self.p_kernel.beta(p_details, p_values, cutoff, magnetic) 
     260        else: 
     261            p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic) 
     262        s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    260263        #print("p_npars",p_npars,s_npars,p_er,s_vr,values[2+p_npars+1:2+p_npars+s_npars]) 
    261264        #call_details.show(values) 
     
    265268        #s_details.show(s_values) 
    266269        #print("=>", s_result) 
    267  
    268         # remember the parts for plotting later 
    269         self.results = [p_result, s_result] 
    270  
    271270        #import pylab as plt 
    272271        #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') 
    273272        #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') 
    274273        #plt.figure() 
    275  
    276         return values[0]*(p_result*s_result) + values[1] 
     274        if beta_mode:#beta 
     275            beta_factor = F1**2/F2 
     276            Sq_eff = 1+beta_factor*(s_result - 1) 
     277            self.results = [F2, Sq_eff,F1,s_result] 
     278            final_result = volfrac*values[0]*(F2 + (F1**2)*(s_result - 1))/volume_avg+values[1] 
     279            #final_result =  volfrac * values[0] * F2 * Sq_eff / volume_avg + values[1] 
     280        else: 
     281            # remember the parts for plotting later 
     282            self.results = [p_result, s_result] 
     283            final_result = values[0]*(p_result*s_result) + values[1] 
     284        return final_result 
    277285 
    278286    def release(self): 
Note: See TracChangeset for help on using the changeset viewer.