Changeset e819f71 in sasmodels


Ignore:
Timestamp:
Mar 3, 2015 3:22:20 PM (10 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
d6adfbe
Parents:
95e861b (diff), 1353f60 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' of github.com:sasview/sasmodels

Files:
5 edited

Legend:

Unmodified
Added
Removed

  • extra/pylint.rc

    r6c8db9e r53d0e24  
    9898 
    9999# Good variable names which should always be accepted, separated by a comma 
    100 good-names=i,j,k,ex,Run,_,x,y,z,qx,qy,qz,n,q,dx,dy,dz,id 
     100good-names=i,j,k,ex,Run,_,x,y,z,qx,qy,qz,n,q,dx,dy,dz,id,Iq,dIq,Qx,Qy,Qz 
    101101 
    102102# Bad variable names which should always be refused, separated by a comma 

  • sasmodels/bumps_model.py

    r9890053 r1353f60  
    22Sasmodels core. 
    33""" 
    4 import sys, os 
    54import datetime 
    65 
    … …  
    98# CRUFT python 2.6 
    109if not hasattr(datetime.timedelta, 'total_seconds'): 
    11     def delay(dt): return dt.days*86400 + dt.seconds + 1e-6*dt.microseconds 
     10    def delay(dt): return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds 
    1211else: 
    1312    def delay(dt): return dt.total_seconds() 
    … …  
    1716try: 
    1817    from .kernelcl import load_model as _loader 
    19 except RuntimeError,exc: 
     18except RuntimeError, exc: 
    2019    import warnings 
    2120    warnings.warn(str(exc)) 
    … …  
    2726    Load model by name. 
    2827    """ 
    29     sasmodels = __import__('sasmodels.models.'+modelname) 
     28    sasmodels = __import__('sasmodels.models.' + modelname) 
    3029    module = getattr(sasmodels.models, modelname, None) 
    3130    model = _loader(module, dtype=dtype) 
    … …  
    4140    """ 
    4241    then = datetime.datetime.now() 
    43     return lambda: delay(datetime.datetime.now()-then) 
     42    return lambda: delay(datetime.datetime.now() - then) 
    4443 
    4544 
    … …  
    5251    data = loader.load(filename) 
    5352    if data is None: 
    54         raise IOError("Data %r could not be loaded"%filename) 
     53        raise IOError("Data %r could not be loaded" % filename) 
    5554    return data 
    5655 
    … …  
    6564    from sas.dataloader.data_info import Data1D 
    6665 
    67     Iq = 100*np.ones_like(q) 
     66    Iq = 100 * np.ones_like(q) 
    6867    dIq = np.sqrt(Iq) 
    69     data = Data1D(q, Iq, dx=0.05*q, dy=dIq) 
     68    data = Data1D(q, Iq, dx=0.05 * q, dy=dIq) 
    7069    data.filename = "fake data" 
    7170    data.qmin, data.qmax = q.min(), q.max() 
    … …  
    8584    if qy is None: 
    8685        qy = qx 
    87     Qx,Qy = np.meshgrid(qx,qy) 
    88     Qx,Qy = Qx.flatten(), Qy.flatten() 
    89     Iq = 100*np.ones_like(Qx) 
     86    Qx, Qy = np.meshgrid(qx, qy) 
     87    Qx, Qy = Qx.flatten(), Qy.flatten() 
     88    Iq = 100 * np.ones_like(Qx) 
    9089    dIq = np.sqrt(Iq) 
    9190    mask = np.ones(len(Iq), dtype='bool') 
    … …  
    10099 
    101100    # 5% dQ/Q resolution 
    102     data.dqx_data = 0.05*Qx 
    103     data.dqy_data = 0.05*Qy 
     101    data.dqx_data = 0.05 * Qx 
     102    data.dqy_data = 0.05 * Qy 
    104103 
    105104    detector = Detector() 
    … …  
    114113    data.Q_unit = "1/A" 
    115114    data.I_unit = "1/cm" 
    116     data.q_data = np.sqrt(Qx**2 + Qy**2) 
     115    data.q_data = np.sqrt(Qx ** 2 + Qy ** 2) 
    117116    data.xaxis("Q_x", "A^{-1}") 
    118117    data.yaxis("Q_y", "A^{-1}") 
    … …  
    129128        data.mask = Ringcut(0, radius)(data) 
    130129        if outer is not None: 
    131             data.mask += Ringcut(outer,np.inf)(data) 
     130            data.mask += Ringcut(outer, np.inf)(data) 
    132131    else: 
    133         data.mask = (data.x>=radius) 
     132        data.mask = (data.x >= radius) 
    134133        if outer is not None: 
    135             data.mask &= (data.x<outer) 
     134            data.mask &= (data.x < outer) 
    136135 
    137136 
    … …  
    165164    import matplotlib.pyplot as plt 
    166165    if hasattr(data, 'qx_data'): 
    167         iq = iq+0 
     166        iq = iq + 0 
    168167        valid = np.isfinite(iq) 
    169168        if scale == 'log': 
    170169            valid[valid] = (iq[valid] > 0) 
    171170            iq[valid] = np.log10(iq[valid]) 
    172         iq[~valid|data.mask] = 0 
     171        iq[~valid | data.mask] = 0 
    173172        #plottable = iq 
    174         plottable = masked_array(iq, ~valid|data.mask) 
     173        plottable = masked_array(iq, ~valid | data.mask) 
    175174        xmin, xmax = min(data.qx_data), max(data.qx_data) 
    176175        ymin, ymax = min(data.qy_data), max(data.qy_data) 
    177176        try: 
    178             if vmin is None: vmin = iq[valid&~data.mask].min() 
    179             if vmax is None: vmax = iq[valid&~data.mask].max() 
     177            if vmin is None: vmin = iq[valid & ~data.mask].min() 
     178            if vmax is None: vmax = iq[valid & ~data.mask].max() 
    180179        except: 
    181180            vmin, vmax = 0, 1 
    182         plt.imshow(plottable.reshape(128,128), 
     181        plt.imshow(plottable.reshape(128, 128), 
    183182                   interpolation='nearest', aspect=1, origin='upper', 
    184183                   extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
    … …  
    189188        else: 
    190189            idx = np.isfinite(iq) 
    191             idx[idx] = (iq[idx]>0) 
     190            idx[idx] = (iq[idx] > 0) 
    192191            plt.loglog(data.x[idx], iq[idx]) 
    193192 
    … …  
    206205        mdata[mdata <= 0] = masked 
    207206    mtheory = masked_array(theory, mdata.mask) 
    208     mresid = masked_array((theory-data.y)/data.dy, mdata.mask) 
     207    mresid = masked_array((theory - data.y) / data.dy, mdata.mask) 
    209208 
    210209    plt.subplot(121) 
    … …  
    214213    plt.subplot(122) 
    215214    plt.plot(data.x, mresid, 'x') 
    216     #plt.axhline(1, color='black', ls='--',lw=1, hold=True) 
    217     #plt.axhline(0, color='black', lw=1, hold=True) 
    218     #plt.axhline(-1, color='black', ls='--',lw=1, hold=True) 
    219215 
    220216def _plot_sesans(data, theory, view): 
    221217    import matplotlib.pyplot as plt 
    222     resid = (theory - data.y)/data.dy 
     218    resid = (theory - data.y) / data.dy 
    223219    plt.subplot(121) 
    224220    plt.errorbar(data.x, data.y, yerr=data.dy) 
    … …  
    236232    """ 
    237233    import matplotlib.pyplot as plt 
    238     resid = (theory-data.data)/data.err_data 
     234    resid = (theory - data.data) / data.err_data 
    239235    plt.subplot(131) 
    240236    plot_data(data, data.data, scale=view) 
    … …  
    270266        self.model = model 
    271267        self.cutoff = cutoff 
    272 # TODO       if  isinstance(data,SESANSData1D)         
     268# TODO       if  isinstance(data,SESANSData1D) 
    273269        if hasattr(data, 'lam'): 
    274270            self.data_type = 'sesans' 
    … …  
    283279        if self.data_type == 'sesans': 
    284280            q = sesans.make_q(data.sample.zacceptance, data.Rmax) 
    285             self.index = slice(None,None) 
     281            self.index = slice(None, None) 
    286282            self.iq = data.y 
    287283            self.diq = data.dy 
    … …  
    289285            q_vectors = [q] 
    290286        elif self.data_type == 'Iqxy': 
    291             self.index = (data.mask==0) & (~np.isnan(data.data)) 
     287            self.index = (data.mask == 0) & (~np.isnan(data.data)) 
    292288            self.iq = data.data[self.index] 
    293289            self.diq = data.err_data[self.index] 
    294290            self._theory = np.zeros_like(data.data) 
    295291            if not partype['orientation'] and not partype['magnetic']: 
    296                 q_vectors = [np.sqrt(data.qx_data**2+data.qy_data**2)] 
     292                q_vectors = [np.sqrt(data.qx_data ** 2 + data.qy_data ** 2)] 
    297293            else: 
    298294                q_vectors = [data.qx_data, data.qy_data] 
    299295        elif self.data_type == 'Iq': 
    300             self.index = (data.x>=data.qmin) & (data.x<=data.qmax) & ~np.isnan(data.y) 
     296            self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y) 
    301297            self.iq = data.y[self.index] 
    302298            self.diq = data.dy[self.index] 
    … …  
    319315            pars.append(name) 
    320316        for name in partype['pd-2d']: 
    321             for xpart,xdefault,xlimits in [ 
     317            for xpart, xdefault, xlimits in [ 
    322318                    ('_pd', 0, limits), 
    323                     ('_pd_n', 35, (0,1000)), 
     319                    ('_pd_n', 35, (0, 1000)), 
    324320                    ('_pd_nsigma', 3, (0, 10)), 
    325321                    ('_pd_type', 'gaussian', None), 
    326322                ]: 
    327                 xname = name+xpart 
     323                xname = name + xpart 
    328324                xvalue = kw.pop(xname, xdefault) 
    329325                if xlimits is not None: 
    … …  
    333329        self._parameter_names = pars 
    334330        if kw: 
    335             raise TypeError("unexpected parameters: %s"%(", ".join(sorted(kw.keys())))) 
     331            raise TypeError("unexpected parameters: %s" % (", ".join(sorted(kw.keys())))) 
    336332        self.update() 
    337333 
    … …  
    340336 
    341337    def numpoints(self): 
     338        """ 
     339            Return the number of points 
     340        """ 
    342341        return len(self.iq) 
    343342 
    344343    def parameters(self): 
    345         return dict((k,getattr(self,k)) for k in self._parameter_names) 
     344        """ 
     345            Return a dictionary of parameters 
     346        """ 
     347        return dict((k, getattr(self, k)) for k in self._parameter_names) 
    346348 
    347349    def theory(self): 
    348350        if 'theory' not in self._cache: 
    349351            if self._fn is None: 
    350                 input = self.model.make_input(self._fn_inputs) 
    351                 self._fn = self.model(input) 
    352  
    353             fixed_pars = [getattr(self,p).value for p in self._fn.fixed_pars] 
     352                input_value = self.model.make_input(self._fn_inputs) 
     353                self._fn = self.model(input_value) 
     354 
     355            fixed_pars = [getattr(self, p).value for p in self._fn.fixed_pars] 
    354356            pd_pars = [self._get_weights(p) for p in self._fn.pd_pars] 
    355357            #print fixed_pars,pd_pars 
    … …  
    357359            #self._theory[:] = self._fn.eval(pars, pd_pars) 
    358360            if self.data_type == 'sesans': 
    359                 P = sesans.hankel(self.data.x, self.data.lam*1e-9, 
    360                                   self.data.sample.thickness/10, self._fn_inputs[0], 
     361                P = sesans.hankel(self.data.x, self.data.lam * 1e-9, 
     362                                  self.data.sample.thickness / 10, self._fn_inputs[0], 
    361363                                  self._theory) 
    362364                self._cache['theory'] = P 
    … …  
    367369    def residuals(self): 
    368370        #if np.any(self.err ==0): print "zeros in err" 
    369         return (self.theory()[self.index]-self.iq)/self.diq 
     371        return (self.theory()[self.index] - self.iq) / self.diq 
    370372 
    371373    def nllf(self): 
    372374        R = self.residuals() 
    373375        #if np.any(np.isnan(R)): print "NaN in residuals" 
    374         return 0.5*np.sum(R**2) 
     376        return 0.5 * np.sum(R ** 2) 
    375377 
    376378    def __call__(self): 
    377         return 2*self.nllf()/self.dof 
     379        return 2 * self.nllf() / self.dof 
    378380 
    379381    def plot(self, view='log'): 
    … …  
    396398            noise = self.diq[self.index] 
    397399        else: 
    398             noise = 0.01*noise 
     400            noise = 0.01 * noise 
    399401            self.diq[self.index] = noise 
    400402        y = self.theory() 
    401         y += y*np.random.randn(*y.shape)*noise 
     403        y += y * np.random.randn(*y.shape) * noise 
    402404        if self.data_type == 'Iq': 
    403405            self.data.y[self.index] = y 
    … …  
    413415 
    414416    def _get_weights(self, par): 
     417        """ 
     418            Get parameter dispersion weights 
     419        """ 
    415420        from . import weights 
    416421 
    417422        relative = self.model.info['partype']['pd-rel'] 
    418423        limits = self.model.info['limits'] 
    419         disperser,value,npts,width,nsigma = [getattr(self, par+ext) 
    420                 for ext in ('_pd_type','','_pd_n','_pd','_pd_nsigma')] 
    421         v,w = weights.get_weights( 
     424        disperser, value, npts, width, nsigma = \ 
     425            [getattr(self, par + ext) for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')] 
     426        v, w = weights.get_weights( 
    422427            disperser, int(npts.value), width.value, nsigma.value, 
    423428            value.value, limits[par], par in relative) 
    424         return v,w/w.max() 
     429        return v, w / w.max() 
    425430 
    426431    def __getstate__(self): 

  • sasmodels/models/stickyhardsphere.py

    rbfb195e r1353f60  
    11# Note: model title and parameter table are inserted automatically 
    2 r"""This calculates the interparticle structure factor for a hard sphere fluid with a narrow attractive well. A perturbative 
    3 solution of the Percus-Yevick closure is used. The strength of the attractive well is described in terms of "stickiness" 
     2r"""This calculates the interparticle structure factor for a hard sphere fluid with 
     3a narrow attractive well. A perturbative solution of the Percus-Yevick closure is used. 
     4The strength of the attractive well is described in terms of "stickiness" 
    45as defined below. The returned value is a dimensionless structure factor, *S(q)*. 
    56 
    6 The perturb (perturbation parameter), |epsilon|, should be held between 0.01 and 0.1. It is best to hold the 
    7 perturbation parameter fixed and let the "stickiness" vary to adjust the interaction strength. The stickiness, |tau|, 
    8 is defined in the equation below and is a function of both the perturbation parameter and the interaction strength. 
    9 |tau| and |epsilon| are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ), the width of the square 
    10 well, |bigdelta| (same units as *R*), and the depth of the well, *Uo*, in units of kT. From the definition, it is clear 
    11 that smaller |tau| means stronger attraction. 
     7The perturb (perturbation parameter), |epsilon|, should be held between 0.01 and 0.1. 
     8It is best to hold the perturbation parameter fixed and let the "stickiness" vary to 
     9adjust the interaction strength. The stickiness, |tau|, is defined in the equation 
     10below and is a function of both the perturbation parameter and the interaction strength. 
     11|tau| and |epsilon| are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ), 
     12the width of the square well, |bigdelta| (same units as *R*), and the depth of the well, 
     13*Uo*, in units of kT. From the definition, it is clear that smaller |tau| means stronger 
     14attraction. 
    1215 
    1316.. image:: img/stickyhardsphere_228.PNG 
    … …  
    1720.. image:: img/stickyhardsphere_229.PNG 
    1821 
    19 The Percus-Yevick (PY) closure was used for this calculation, and is an adequate closure for an attractive interparticle 
    20 potential. This solution has been compared to Monte Carlo simulations for a square well fluid, with good agreement. 
     22The Percus-Yevick (PY) closure was used for this calculation, and is an adequate closure 
     23for an attractive interparticle potential. This solution has been compared to Monte Carlo 
     24simulations for a square well fluid, with good agreement. 
    2125 
    22 The true particle volume fraction, |phi|, is not equal to *h*, which appears in most of the reference. The two are 
    23 related in equation (24) of the reference. The reference also describes the relationship between this perturbation 
    24 solution and the original sticky hard sphere (or adhesive sphere) model by Baxter. 
     26The true particle volume fraction, |phi|, is not equal to *h*, which appears in most of 
     27the reference. The two are related in equation (24) of the reference. The reference also 
     28describes the relationship between this perturbation solution and the original sticky hard 
     29sphere (or adhesive sphere) model by Baxter. 
    2530 
    26 NB: The calculation can go haywire for certain combinations of the input parameters, producing unphysical solutions - in 
    27 this case errors are reported to the command window and the *S(q)* is set to -1 (so it will disappear on a log-log 
    28 plot). Use tight bounds to keep the parameters to values that you know are physical (test them) and keep nudging them 
    29 until the optimization does not hit the constraints. 
     31NB: The calculation can go haywire for certain combinations of the input parameters, 
     32producing unphysical solutions - in this case errors are reported to the command window and 
     33the *S(q)* is set to -1 (so it will disappear on a log-log plot). Use tight bounds to keep 
     34the parameters to values that you know are physical (test them) and keep nudging them until 
     35the optimization does not hit the constraints. 
    3036 
    31 In sasview the effective radius will be calculated from the parameters used in the form factor P(Q) that this  
    32 S(Q) is combined with. 
     37In sasview the effective radius will be calculated from the parameters used in the 
     38form factor P(Q) that this S(Q) is combined with. 
    3339 
    34 For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as 
     40For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where 
     41the *q* vector is defined as 
    3542 
    3643.. math:: 
    … …  
    5764 
    5865# TODO: refactor so that we pull in the old sansmodels.c_extensions 
    59   
    60 from numpy import pi, inf 
     66 
     67from numpy import inf 
    6168 
    6269name = "stickyhardsphere" 
    … …  
    6471description = """\ 
    6572        [Sticky hard sphere structure factor, with Percus-Yevick closure] 
    66         Interparticle structure factor S(Q)for a hard sphere fluid with  
     73        Interparticle structure factor S(Q)for a hard sphere fluid with 
    6774                a narrow attractive well. Fits are prone to deliver non-physical 
    68                 parameters, use with care and read the references in the full manual.  
    69                 In sasview the effective radius will be calculated from the  
     75                parameters, use with care and read the references in the full manual. 
     76                In sasview the effective radius will be calculated from the 
    7077                parameters used in P(Q). 
    7178""" 
    … …  
    7380 
    7481parameters = [ 
    75 #   [ "name", "units", default, [lower, upper], "type", 
    76 #     "description" ], 
    77     [ "effect_radius", "Ang", 50.0, [0, inf], "volume", 
    78       "effective radius of hard sphere" ], 
    79     [ "volfraction", "", 0.2, [0, 0.74], "", 
    80       "volume fraction of hard spheres" ], 
    81     [ "perturb", "", 0.05, [0.01, 0.1], "", 
    82       "perturbation parameter, epsilon" ], 
    83     [ "stickiness", "",  0.20, [-inf,inf], "", 
    84       "stickiness, tau" ], 
     82    #   [ "name", "units", default, [lower, upper], "type", 
     83    #     "description" ], 
     84    ["effect_radius", "Ang", 50.0, [0, inf], "volume", 
     85     "effective radius of hard sphere"], 
     86    ["volfraction", "", 0.2, [0, 0.74], "", 
     87     "volume fraction of hard spheres"], 
     88    ["perturb", "", 0.05, [0.01, 0.1], "", 
     89     "perturbation parameter, epsilon"], 
     90    ["stickiness", "", 0.20, [-inf, inf], "", 
     91     "stickiness, tau"], 
    8592    ] 
    86          
     93 
    8794# No volume normalization despite having a volume parameter 
    8895# This should perhaps be volume normalized? 
    … …  
    96103        double lam,lam2,test,mu,alpha,beta; 
    97104        double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq; 
    98          
     105 
    99106        onemineps = 1.0-perturb; 
    100107        eta = volfraction/onemineps/onemineps/onemineps; 
    101          
     108 
    102109        sig = 2.0 * effect_radius; 
    103110        aa = sig/onemineps; 
    … …  
    153160        // 
    154161        sq = 1.0/(aq*aq +bq*bq); 
    155          
     162 
    156163        return(sq); 
    157164""" 
    … …  
    166173oldname = 'StickyHSStructure' 
    167174oldpars = dict() 
    168 demo = dict(effect_radius = 200,volfraction = 0.2,perturb=0.05,stickiness=0.2,effect_radius_pd = 0.1,effect_radius_pd_n = 40) 
     175demo = dict(effect_radius=200, volfraction=0.2, perturb=0.05, 
     176            stickiness=0.2, effect_radius_pd=0.1, effect_radius_pd_n=40) 
    169177 
    170178 

  • sasmodels/kernel_template.c

    rf734e7d r95e861b  
    1010#ifndef USE_OPENCL 
    1111#  ifdef __cplusplus 
     12     #include <cstdio> 
    1213     #include <cmath> 
    1314     #if defined(_MSC_VER) 
    … …  
    2021       { svar=sin(angle); cvar=cos(angle); } 
    2122#  else 
     23     #include <stdio.h> 
    2224     #include <math.h> 
    2325     #if defined(_MSC_VER) 

  • sasmodels/models/lamellarCailleHG_kernel.c

    rcd55ac3 r95e861b  
    5555    //temp *= cos(dd*qval*ii/(1.0+t1)); 
    5656    temp *= cos(dd*qval*ii); 
     57    //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii); 
    5758    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 
    5859    temp *= exp(-t2/2.0 ); 
    … …  
    6869 
    6970  inten *= 1.0e-04;   // 1/A to 1/cm 
    70  
    7171  return(inten); 
    7272} 
Note: See TracChangeset for help on using the changeset viewer.