Changeset 9acade6 in sasmodels


Ignore:
Timestamp:
Sep 27, 2016 7:52:14 AM (8 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:
bb46723
Parents:
6cefbc9 (diff), 5fd684d (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

Location:
sasmodels
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/models/pringle.c

    rc047acf r58c3367  
    4242    } 
    4343 
    44     #if 1 
    45     // TODO: should the normalization be to R^2 or the last value, r^2 
    46     // the gaussian window does not go all the way from 0 to 1. 
    47     //radius = Gauss76Z[75] * zm + zb; 
    48     *Sn = zm*sumS / (r*r); 
    49     *Cn = zm*sumC / (r*r); 
    50     #else 
    5144    *Sn = zm*sumS / (radius*radius); 
    5245    *Cn = zm*sumC / (radius*radius); 
    53     #endif 
    5446} 
    5547 
     
    109101double form_volume(double radius, double thickness, double alpha, double beta) 
    110102{ 
    111     // TODO: Normalize by form volume 
    112     //return M_PI*radius*radius*thickness; 
    113     return 1.0; 
     103    return M_PI*radius*radius*thickness; 
    114104} 
    115105 
     
    126116    double contrast = sld_pringle - sld_solvent; 
    127117    double volume = M_PI*radius*radius*thickness; 
    128     // TODO: If normalize by form volume, need an extra volume here 
    129     //return 1.0e-4*form * square(contrast * volume); 
    130     return 1.0e-4*form * square(contrast) * volume; 
     118    return 1.0e-4*form * square(contrast * volume); 
    131119} 
  • sasmodels/models/pringle.py

    r40a87fa r5fd684d  
    1717.. math:: 
    1818 
    19     C_n = \int^{R}_{0} r dr\cos(qr^2\alpha \cos{\psi}) 
     19    C_n = \frac{1}{r^2}\int^{R}_{0} r dr\cos(qr^2\alpha \cos{\psi}) 
    2020    J_n\left( qr^2\beta \cos{\psi}\right) 
    2121    J_{2n}\left( qr \sin{\psi}\right) 
     
    2323.. math:: 
    2424 
    25     S_n = \int^{R}_{0} r dr\sin(qr^2\alpha \cos{\psi}) 
     25    S_n = \frac{1}{r^2}\int^{R}_{0} r dr\sin(qr^2\alpha \cos{\psi}) 
    2626    J_n\left( qr^2\beta \cos{\psi}\right) 
    2727    J_{2n}\left( qr \sin{\psi}\right) 
     
    4747**Last Modified by:** Wojciech Wpotrzebowski **on:** March 20, 2016 
    4848 
    49 **Last Reviewed by:** Paul Butler **on:** March 21, 2016 
     49**Last Reviewed by:** Andrew Jackson **on:** September 26, 2016 
    5050 
    5151""" 
     
    102102      'sld_pringle': 1.0, 
    103103      'sld_solvent': 6.3, 
    104       'background': 6.3, 
    105      }, 0.1, 16.185532], 
     104      'background': 0.001, 
     105     }, 0.1, 9.87676], 
    106106 
    107107    [{'scale' : 1.0, 
     
    112112      'sld_pringle': 1.0, 
    113113      'sld_solvent': 6.3, 
    114       'background': 6.3, 
    115      }, 0.01, 297.153496], 
     114      'background': 0.001, 
     115     }, 0.01, 290.56723], 
    116116 
    117117    [{'scale' : 1.0, 
     
    122122      'sld_pringle': 1.0, 
    123123      'sld_solvent': 6.3, 
    124       'background': 6.3, 
    125      }, 0.001, 324.021256415], 
    126  
    127     [{'scale' : 1.0, 
    128       'radius': 60.0, 
    129       'thickness': 10.0, 
    130       'alpha': 0.001, 
    131       'beta': 0.02, 
    132       'sld_pringle': 1.0, 
    133       'sld_solvent': 6.3, 
    134       'background': 6.3, 
    135      }, (0.001, 90.0), 6.30000026876], 
     124      'background': 0.001, 
     125     }, 0.001, 317.40847], 
    136126] 
  • sasmodels/sasview_model.py

    r3bcb88c r9c1a59c  
    565565        parameters = self._model_info.parameters 
    566566        pairs = [self._get_weights(p) for p in parameters.call_parameters] 
     567        #weights.plot_weights(self._model_info, pairs) 
    567568        call_details, values, is_magnetic = make_kernel_args(calculator, pairs) 
    568569        #call_details.show() 
     
    618619            # remember them is kind of funky. 
    619620            # Note: can't seem to get disperser parameters from sasview 
    620             # (1) Could create a sasview model that has not yet # been 
     621            # (1) Could create a sasview model that has not yet been 
    621622            # converted, assign the disperser to one of its polydisperse 
    622623            # parameters, then retrieve the disperser parameters from the 
    623             # sasview model.  (2) Could write a disperser parameter retriever 
    624             # in sasview.  (3) Could modify sasview to use sasmodels.weights 
    625             # dispersers. 
     624            # sasview model. 
     625            # (2) Could write a disperser parameter retriever in sasview. 
     626            # (3) Could modify sasview to use sasmodels.weights dispersers. 
    626627            # For now, rely on the fact that the sasview only ever uses 
    627628            # new dispersers in the set_dispersion call and create a new 
    628629            # one instead of trying to assign parameters. 
    629             dispersion = weights.MODELS[dispersion.type]() 
    630630            self.dispersion[parameter] = dispersion.get_pars() 
    631631        else: 
     
    658658        elif par.polydisperse: 
    659659            dis = self.dispersion[par.name] 
    660             value, weight = weights.get_weights( 
    661                 dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
    662                 self.params[par.name], par.limits, par.relative_pd) 
     660            if dis['type'] == 'array': 
     661                value, weight = dis['values'], dis['weights'] 
     662            else: 
     663                value, weight = weights.get_weights( 
     664                    dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
     665                    self.params[par.name], par.limits, par.relative_pd) 
    663666            return value, weight / np.sum(weight) 
    664667        else: 
  • sasmodels/weights.py

    r56b2687 r6cefbc9  
    33""" 
    44# TODO: include dispersion docs with the disperser models 
     5from __future__ import division, print_function 
     6 
    57from math import sqrt  # type: ignore 
     8from collections import OrderedDict 
     9 
    610import numpy as np  # type: ignore 
    711from scipy.special import gammaln  # type: ignore 
     
    5458            else: 
    5559                return np.array([], 'd'), np.array([], 'd') 
    56         return self._weights(center, sigma, lb, ub) 
     60        x, px = self._weights(center, sigma, lb, ub) 
     61        return x, px 
    5762 
    5863    def _weights(self, center, sigma, lb, ub): 
     
    7984    default = dict(npts=35, width=0, nsigmas=3) 
    8085    def _weights(self, center, sigma, lb, ub): 
     86        # TODO: sample high probability regions more densely 
     87        # i.e., step uniformly in cumulative density rather than x value 
     88        # so weight = 1/Npts for all weights, but values are unevenly spaced 
    8189        x = self._linspace(center, sigma, lb, ub) 
    8290        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma)) 
     
    165173 
    166174    def _weights(self, center, sigma, lb, ub): 
    167         # TODO: interpolate into the array dispersion using npts 
    168         x = center + self.values*sigma 
     175        # TODO: rebin the array dispersion using npts 
     176        # TODO: use a distribution that can be recentered and scaled 
     177        x = self.values 
     178        #x = center + self.values*sigma 
    169179        idx = (x >= lb) & (x <= ub) 
    170180        x = x[idx] 
     
    174184 
    175185# dispersion name -> disperser lookup table. 
    176 MODELS = dict((d.type, d) for d in ( 
    177     GaussianDispersion, RectangleDispersion, 
    178     ArrayDispersion, SchulzDispersion, LogNormalDispersion 
     186# Maintain order since this is used by sasview GUI to order the options in 
     187# the dispersion type combobox. 
     188MODELS = OrderedDict((d.type, d) for d in ( 
     189    RectangleDispersion, 
     190    ArrayDispersion, 
     191    LogNormalDispersion, 
     192    GaussianDispersion, 
     193    SchulzDispersion, 
    179194)) 
    180195 
     
    194209    *value* is the value of the parameter in the model. 
    195210 
    196     *limits* is [lb, ub], the lower and upper bound of the weight value. 
     211    *limits* is [lb, ub], the lower and upper bound on the possible values. 
    197212 
    198213    *relative* is true if *width* is defined in proportion to the value 
     
    201216    Returns *(value, weight)*, where *value* and *weight* are vectors. 
    202217    """ 
     218    if disperser == "array": 
     219        raise NotImplementedError("Don't handle arrays through get_weights; use values and weights directly") 
    203220    cls = MODELS[disperser] 
    204221    obj = cls(n, width, nsigmas) 
    205222    v, w = obj.get_weights(value, limits[0], limits[1], relative) 
    206223    return v, w 
     224 
     225 
     226def plot_weights(model_info, pairs): 
     227    # type: (ModelInfo, List[Tuple[np.ndarray, np.ndarray]]) -> None 
     228    """ 
     229    Plot the weights returned by :func:`get_weights`. 
     230 
     231    *model_info* is 
     232    :param model_info: 
     233    :param pairs: 
     234    :return: 
     235    """ 
     236    import pylab 
     237 
     238    if any(len(values)>1 for values, weights in pairs): 
     239        labels = [p.name for p in model_info.parameters.call_parameters] 
     240        pylab.interactive(True) 
     241        pylab.figure() 
     242        for (v,w), s in zip(pairs, labels): 
     243            if len(v) > 1: 
     244                #print("weights for", s, v, w) 
     245                pylab.plot(v, w, '-o', label=s) 
     246        pylab.grid(True) 
     247        pylab.legend() 
Note: See TracChangeset for help on using the changeset viewer.