Changeset 13bd583 in sasmodels


Ignore:
Timestamp:
Feb 12, 2018 8:24:15 PM (7 years ago)
Author:
GitHub <noreply@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
a2ca6e5
Parents:
e309e23 (diff), f6fd413 (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.
git-author:
Paul Butler <butlerpd@…> (02/12/18 20:24:15)
git-committer:
GitHub <noreply@…> (02/12/18 20:24:15)
Message:

Merge pull request #62 from SasView?/pytest

add support for pytest and use it on travis/appveyor

Files:
5 added
14 edited

Legend:

Unmodified
Added
Removed

  • .travis.yml

    r2d09df1 r335271e  
    3737- conda update --yes conda 
    3838- conda info -a 
    39 - conda install --yes python=$PY numpy scipy cython mako cffi 
     39- conda install --yes python=$PY numpy scipy matplotlib docutils setuptools pytest 
    4040install: 
    4141- pip install bumps 
    42 - pip install unittest-xml-reporting 
     42- pip install unittest-xml-reporting tinycc 
    4343script: 
    4444- python --version 
    45 - python -m sasmodels.model_test -v dll all 
     45#- python -m sasmodels.model_test -v dll all 
     46#- python -m pytest -v --cache-clear 
     47- python setup.py test --pytest-args -v 
    4648before_deploy: 
    4749- echo -e "Host danse.chem.utk.edu\n\tStrictHostKeyChecking no\n" >> ~/.ssh/config 

  • appveyor.yml

    rfefc185 r335271e  
    4444 
    4545install: 
    46     # Set the CONDA_NPY, although it has no impact on the actual build.  
     46    # Set the CONDA_NPY, although it has no impact on the actual build. 
    4747    # We need this because of a test within conda-build. 
    4848    - cmd: set CONDA_NPY=19 
    … …  
    6868    #- cmd: conda install --yes --quiet obvious-ci 
    6969    # 2018-01-19 PAK: skipping toolchain cython and cffi (these were for pyopencl?) 
    70     - cmd: conda install --yes --quiet numpy scipy 
     70    - cmd: conda install --yes --quiet numpy scipy matplotlib docutils setuptools pytest 
    7171    # 2018-01-19 PAK: skipping pyopencl; this would be needed for deploy but maybe not for test 
    7272    #- cmd: conda install --yes --channel conda-forge pyopencl 
    … …  
    8686    #- "%CMD_IN_ENV% python -m sasmodels.model_test dll all" 
    8787    #- cmd /E:ON /V:ON /C obvci_appveyor_python_build_env.cmd python -m sasmoels.model_test dll all 
    88     - python -m sasmodels.model_test dll all 
     88    #- python -m sasmodels.model_test dll all 
     89    #- nosetests -v sasmodels/*.py 
     90    #- python -m pytest -v --cache-clear 
     91    - python setup.py test --pytest-args -v 

  • sasmodels/compare.py

    r2a7e20e r3221de0  
    4040from . import core 
    4141from . import kerneldll 
     42from . import kernelcl 
    4243from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4344from .direct_model import DirectModel, get_mesh 
    … …  
    623624 
    624625 
     626def time_calculation(calculator, pars, evals=1): 
     627    # type: (Calculator, ParameterSet, int) -> Tuple[np.ndarray, float] 
     628    """ 
     629    Compute the average calculation time over N evaluations. 
     630 
     631    An additional call is generated without polydispersity in order to 
     632    initialize the calculation engine, and make the average more stable. 
     633    """ 
     634    # initialize the code so time is more accurate 
     635    if evals > 1: 
     636        calculator(**suppress_pd(pars)) 
     637    toc = tic() 
     638    # make sure there is at least one eval 
     639    value = calculator(**pars) 
     640    for _ in range(evals-1): 
     641        value = calculator(**pars) 
     642    average_time = toc()*1000. / evals 
     643    #print("I(q)",value) 
     644    return value, average_time 
     645 
     646def make_data(opts): 
     647    # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray] 
     648    """ 
     649    Generate an empty dataset, used with the model to set Q points 
     650    and resolution. 
     651 
     652    *opts* contains the options, with 'qmax', 'nq', 'res', 
     653    'accuracy', 'is2d' and 'view' parsed from the command line. 
     654    """ 
     655    qmin, qmax, nq, res = opts['qmin'], opts['qmax'], opts['nq'], opts['res'] 
     656    if opts['is2d']: 
     657        q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray 
     658        data = empty_data2D(q, resolution=res) 
     659        data.accuracy = opts['accuracy'] 
     660        set_beam_stop(data, qmin) 
     661        index = ~data.mask 
     662    else: 
     663        if opts['view'] == 'log' and not opts['zero']: 
     664            q = np.logspace(math.log10(qmin), math.log10(qmax), nq) 
     665        else: 
     666            q = np.linspace(qmin, qmax, nq) 
     667        if opts['zero']: 
     668            q = np.hstack((0, q)) 
     669        data = empty_data1D(q, resolution=res) 
     670        index = slice(None, None) 
     671    return data, index 
     672 
    625673DTYPE_MAP = { 
    626674    'half': '16', 
    … …  
    643691    Return a model calculator using the OpenCL calculation engine. 
    644692    """ 
    645     if not core.HAVE_OPENCL: 
    646         raise RuntimeError("OpenCL not available") 
    647     model = core.build_model(model_info, dtype=dtype, platform="ocl") 
    648     calculator = DirectModel(data, model, cutoff=cutoff) 
    649     calculator.engine = "OCL%s"%DTYPE_MAP[str(model.dtype)] 
    650     return calculator 
    651693 
    652694def eval_ctypes(model_info, data, dtype='double', cutoff=0.): 
    … …  
    660702    return calculator 
    661703 
    662 def time_calculation(calculator, pars, evals=1): 
    663     # type: (Calculator, ParameterSet, int) -> Tuple[np.ndarray, float] 
    664     """ 
    665     Compute the average calculation time over N evaluations. 
    666  
    667     An additional call is generated without polydispersity in order to 
    668     initialize the calculation engine, and make the average more stable. 
    669     """ 
    670     # initialize the code so time is more accurate 
    671     if evals > 1: 
    672         calculator(**suppress_pd(pars)) 
    673     toc = tic() 
    674     # make sure there is at least one eval 
    675     value = calculator(**pars) 
    676     for _ in range(evals-1): 
    677         value = calculator(**pars) 
    678     average_time = toc()*1000. / evals 
    679     #print("I(q)",value) 
    680     return value, average_time 
    681  
    682 def make_data(opts): 
    683     # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray] 
    684     """ 
    685     Generate an empty dataset, used with the model to set Q points 
    686     and resolution. 
    687  
    688     *opts* contains the options, with 'qmax', 'nq', 'res', 
    689     'accuracy', 'is2d' and 'view' parsed from the command line. 
    690     """ 
    691     qmin, qmax, nq, res = opts['qmin'], opts['qmax'], opts['nq'], opts['res'] 
    692     if opts['is2d']: 
    693         q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray 
    694         data = empty_data2D(q, resolution=res) 
    695         data.accuracy = opts['accuracy'] 
    696         set_beam_stop(data, qmin) 
    697         index = ~data.mask 
    698     else: 
    699         if opts['view'] == 'log' and not opts['zero']: 
    700             q = np.logspace(math.log10(qmin), math.log10(qmax), nq) 
    701         else: 
    702             q = np.linspace(qmin, qmax, nq) 
    703         if opts['zero']: 
    704             q = np.hstack((0, q)) 
    705         data = empty_data1D(q, resolution=res) 
    706         index = slice(None, None) 
    707     return data, index 
    708  
    709704def make_engine(model_info, data, dtype, cutoff, ngauss=0): 
    710705    # type: (ModelInfo, Data, str, float) -> Calculator 
    … …  
    718713        set_integration_size(model_info, ngauss) 
    719714 
    720     if dtype is None or not dtype.endswith('!'): 
    721         return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) 
    722     else: 
    723         return eval_ctypes(model_info, data, dtype=dtype[:-1], cutoff=cutoff) 
     715    if dtype != "default" and not dtype.endswith('!') and not kernelcl.use_opencl(): 
     716        raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR) 
     717 
     718    model = core.build_model(model_info, dtype=dtype, platform="ocl") 
     719    calculator = DirectModel(data, model, cutoff=cutoff) 
     720    engine_type = calculator._model.__class__.__name__.replace('Model','').upper() 
     721    bits = calculator._model.dtype.itemsize*8 
     722    precision = "fast" if getattr(calculator._model, 'fast', False) else str(bits) 
     723    calculator.engine = "%s[%s]" % (engine_type, precision) 
     724    return calculator 
    724725 
    725726def _show_invalid(data, theory): 

  • sasmodels/convert.py

    r2d81cfe ra69d8cd  
    633633def test_backward_forward(): 
    634634    from .core import list_models 
     635    L = lambda name: _check_one(name, seed=1) 
    635636    for name in list_models('all'): 
    636         L = lambda: _check_one(name, seed=1) 
    637         L.description = name 
    638         yield L 
     637        yield L, name 

  • sasmodels/core.py

    r7a516d0 r3221de0  
    2121from . import mixture 
    2222from . import kernelpy 
     23from . import kernelcl 
    2324from . import kerneldll 
    2425from . import custom 
    25  
    26 if os.environ.get("SAS_OPENCL", "").lower() == "none": 
    27     HAVE_OPENCL = False 
    28 else: 
    29     try: 
    30         from . import kernelcl 
    31         HAVE_OPENCL = True 
    32     except Exception: 
    33         HAVE_OPENCL = False 
    34  
    35 CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "") 
    36 if CUSTOM_MODEL_PATH == "": 
    37     CUSTOM_MODEL_PATH = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models") 
    38     if not os.path.isdir(CUSTOM_MODEL_PATH): 
    39         os.makedirs(CUSTOM_MODEL_PATH) 
    4026 
    4127# pylint: disable=unused-import 
    … …  
    4733    pass 
    4834# pylint: enable=unused-import 
     35 
     36CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "") 
     37if CUSTOM_MODEL_PATH == "": 
     38    CUSTOM_MODEL_PATH = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models") 
     39    if not os.path.isdir(CUSTOM_MODEL_PATH): 
     40        os.makedirs(CUSTOM_MODEL_PATH) 
    4941 
    5042# TODO: refactor composite model support 
    … …  
    276268    if platform is None: 
    277269        platform = "ocl" 
    278     if platform == "ocl" and not HAVE_OPENCL or not model_info.opencl: 
     270    if not kernelcl.use_opencl() or not model_info.opencl: 
    279271        platform = "dll" 
    280272 
    … …  
    320312    print("\n".join(list_models(kind))) 
    321313 
    322 def test_load(): 
    323     # type: () -> None 
    324     """Check that model load works""" 
     314def test_composite_order(): 
    325315    def test_models(fst, snd): 
    326316        """Confirm that two models produce the same parameters""" 
    327317        fst = load_model(fst) 
    328318        snd = load_model(snd) 
    329         # Remove the upper case characters frin the parameters, since 
    330         # they lasndel the order of events and we specfically are 
    331         # changin that aspect 
     319        # Un-disambiguate parameter names so that we can check if the same 
     320        # parameters are in a pair of composite models. Since each parameter in 
     321        # the mixture model is tagged as e.g., A_sld, we ought to use a 
     322        # regex subsitution s/^[A-Z]+_/_/, but removing all uppercase letters 
     323        # is good enough. 
    332324        fst = [[x for x in p.name if x == x.lower()] for p in fst.info.parameters.kernel_parameters] 
    333325        snd = [[x for x in p.name if x == x.lower()] for p in snd.info.parameters.kernel_parameters] 
    334326        assert sorted(fst) == sorted(snd), "{} != {}".format(fst, snd) 
    335327 
    336  
    337     test_models( 
     328    def build_test(first, second): 
     329        test = lambda description: test_models(first, second) 
     330        description = first + " vs. " + second 
     331        return test, description 
     332 
     333    yield build_test( 
    338334        "cylinder+sphere", 
    339335        "sphere+cylinder") 
    340     test_models( 
     336    yield build_test( 
    341337        "cylinder*sphere", 
    342338        "sphere*cylinder") 
    343     test_models( 
     339    yield build_test( 
    344340        "cylinder@hardsphere*sphere", 
    345341        "sphere*cylinder@hardsphere") 
    346     test_models( 
     342    yield build_test( 
    347343        "barbell+sphere*cylinder@hardsphere", 
    348344        "sphere*cylinder@hardsphere+barbell") 
    349     test_models( 
     345    yield build_test( 
    350346        "barbell+cylinder@hardsphere*sphere", 
    351347        "cylinder@hardsphere*sphere+barbell") 
    352     test_models( 
     348    yield build_test( 
    353349        "barbell+sphere*cylinder@hardsphere", 
    354350        "barbell+cylinder@hardsphere*sphere") 
    355     test_models( 
     351    yield build_test( 
    356352        "sphere*cylinder@hardsphere+barbell", 
    357353        "cylinder@hardsphere*sphere+barbell") 
    358     test_models( 
     354    yield build_test( 
    359355        "barbell+sphere*cylinder@hardsphere", 
    360356        "cylinder@hardsphere*sphere+barbell") 
    361     test_models( 
     357    yield build_test( 
    362358        "barbell+cylinder@hardsphere*sphere", 
    363359        "sphere*cylinder@hardsphere+barbell") 
    364360 
     361def test_composite(): 
     362    # type: () -> None 
     363    """Check that model load works""" 
    365364    #Test the the model produces the parameters that we would expect 
    366365    model = load_model("cylinder@hardsphere*sphere") 

  • sasmodels/jitter.py

    r0d5a655 r42158d2  
    1313import argparse 
    1414 
    15 import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
     15try: # CRUFT: travis-ci does not support mpl_toolkits.mplot3d 
     16    import mpl_toolkits.mplot3d  # Adds projection='3d' option to subplot 
     17except ImportError: 
     18    pass 
     19 
    1620import matplotlib.pyplot as plt 
    1721from matplotlib.widgets import Slider, CheckButtons 

  • sasmodels/kernelcl.py

    r2d81cfe rb4272a2  
    5858import numpy as np  # type: ignore 
    5959 
     60 
     61# Attempt to setup opencl. This may fail if the opencl package is not 
     62# installed or if it is installed but there are no devices available. 
    6063try: 
    61     #raise NotImplementedError("OpenCL not yet implemented for new kernel template") 
    6264    import pyopencl as cl  # type: ignore 
     65    from pyopencl import mem_flags as mf 
     66    from pyopencl.characterize import get_fast_inaccurate_build_options 
    6367    # Ask OpenCL for the default context so that we know that one exists 
    6468    cl.create_some_context(interactive=False) 
     69    HAVE_OPENCL = True 
     70    OPENCL_ERROR = "" 
    6571except Exception as exc: 
    66     warnings.warn("OpenCL startup failed with ***" 
    67                   + str(exc) + "***; using C compiler instead") 
    68     raise RuntimeError("OpenCL not available") 
    69  
    70 from pyopencl import mem_flags as mf 
    71 from pyopencl.characterize import get_fast_inaccurate_build_options 
     72    HAVE_OPENCL = False 
     73    OPENCL_ERROR = str(exc) 
    7274 
    7375from . import generate 
    … …  
    102104        cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS] 
    103105 
    104 fix_pyopencl_include() 
    105  
     106if HAVE_OPENCL: 
     107    fix_pyopencl_include() 
    106108 
    107109# The max loops number is limited by the amount of local memory available 
    … …  
    128130""" 
    129131 
     132def use_opencl(): 
     133    return HAVE_OPENCL and os.environ.get("SAS_OPENCL", "").lower() != "none" 
    130134 
    131135ENV = None 
     136def reset_environment(): 
     137    """ 
     138    Call to create a new OpenCL context, such as after a change to SAS_OPENCL. 
     139    """ 
     140    global ENV 
     141    ENV = GpuEnvironment() if use_opencl() else None 
     142 
    132143def environment(): 
    133144    # type: () -> "GpuEnvironment" 
    … …  
    137148    This provides an OpenCL context and one queue per device. 
    138149    """ 
    139     global ENV 
    140150    if ENV is None: 
    141         ENV = GpuEnvironment() 
     151        if not HAVE_OPENCL: 
     152            raise RuntimeError("OpenCL startup failed with ***" 
     153                            + OPENCL_ERROR + "***; using C compiler instead") 
     154        reset_environment() 
     155        if ENV is None: 
     156            raise RuntimeError("SAS_OPENCL=None in environment") 
    142157    return ENV 
    143158 

  • sasmodels/model_test.py

    r2d81cfe r3221de0  
    6262from .exception import annotate_exception 
    6363from .modelinfo import expand_pars 
     64from .kernelcl import use_opencl 
    6465 
    6566# pylint: disable=unused-import 
    … …  
    123124 
    124125        if is_py:  # kernel implemented in python 
    125             test_name = "Model: %s, Kernel: python"%model_name 
     126            test_name = "%s-python"%model_name 
    126127            test_method_name = "test_%s_python" % model_info.id 
    127128            test = ModelTestCase(test_name, model_info, 
    … …  
    134135 
    135136            # test using dll if desired 
    136             if 'dll' in loaders or not core.HAVE_OPENCL: 
    137                 test_name = "Model: %s, Kernel: dll"%model_name 
     137            if 'dll' in loaders or not use_opencl(): 
     138                test_name = "%s-dll"%model_name 
    138139                test_method_name = "test_%s_dll" % model_info.id 
    139140                test = ModelTestCase(test_name, model_info, 
    … …  
    145146 
    146147            # test using opencl if desired and available 
    147             if 'opencl' in loaders and core.HAVE_OPENCL: 
    148                 test_name = "Model: %s, Kernel: OpenCL"%model_name 
     148            if 'opencl' in loaders and use_opencl(): 
     149                test_name = "%s-opencl"%model_name 
    149150                test_method_name = "test_%s_opencl" % model_info.id 
    150151                # Using dtype=None so that the models that are only 
    … …  
    160161 
    161162    return suite 
    162  
    163163 
    164164def _hide_model_case_from_nose(): 
    … …  
    368368 
    369369    # Build a test suite containing just the model 
    370     loaders = ['opencl'] if core.HAVE_OPENCL else ['dll'] 
     370    loaders = ['opencl'] if use_opencl() else ['dll'] 
    371371    models = [model] 
    372372    try: 
    … …  
    425425        verbosity = 1 
    426426    if models and models[0] == 'opencl': 
    427         if not core.HAVE_OPENCL: 
     427        if not use_opencl(): 
    428428            print("opencl is not available") 
    429429            return 1 
    … …  
    435435        models = models[1:] 
    436436    elif models and models[0] == 'opencl_and_dll': 
    437         loaders = ['opencl', 'dll'] if core.HAVE_OPENCL else ['dll'] 
     437        loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
    438438        models = models[1:] 
    439439    else: 
    440         loaders = ['opencl', 'dll'] if core.HAVE_OPENCL else ['dll'] 
     440        loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
    441441    if not models: 
    442442        print("""\ 
    … …  
    467467    Run "nosetests sasmodels" on the command line to invoke it. 
    468468    """ 
    469     loaders = ['opencl', 'dll'] if core.HAVE_OPENCL else ['dll'] 
     469    loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
    470470    tests = make_suite(loaders, ['all']) 
    471     for test_i in tests: 
    472         # In order for nosetest to see the correct test name, need to set 
    473         # the description attribute of the returned function.  Since we 
    474         # can't do this for the returned instance, wrap it in a lambda and 
    475         # set the description on the lambda.  Otherwise we could just do: 
    476         #    yield test_i.run_all 
    477         L = lambda: test_i.run_all() 
    478         L.description = test_i.test_name 
    479         yield L 
     471    def build_test(test): 
     472        # In order for nosetest to show the test name, wrap the test.run_all 
     473        # instance in function that takes the test name as a parameter which 
     474        # will be displayed when the test is run.  Do this as a function so 
     475        # that it properly captures the context for tests that captured and 
     476        # run later.  If done directly in the for loop, then the looping 
     477        # variable test will be shared amongst all the tests, and we will be 
     478        # repeatedly testing vesicle. 
     479 
     480        # Note: in sasview sas.sasgui.perspectives.fitting.gpu_options 
     481        # requires that the test.description field be set. 
     482        wrap = lambda: test.run_all() 
     483        wrap.description = test.test_name 
     484        return wrap 
     485        # The following would work with nosetests and pytest: 
     486        #     return lambda name: test.run_all(), test.test_name 
     487 
     488    for test in tests: 
     489        yield build_test(test) 
    480490 
    481491 

  • sasmodels/sasview_model.py

    r2d81cfe r3221de0  
    665665 
    666666    def _calculate_Iq(self, qx, qy=None): 
    667         #core.HAVE_OPENCL = False 
    668667        if self._model is None: 
    669668            self._model = core.build_model(self._model_info) 
    … …  
    859858    # type: () -> None 
    860859    """ 
    861     Load and run cylinder model from sas.models.CylinderModel 
     860    Load and run cylinder model as sas-models-CylinderModel 
    862861    """ 
    863862    if not SUPPORT_OLD_STYLE_PLUGINS: 

  • setup.py

    r32e3c9b r1f991d6  
    1 try: 
    2     from setuptools import setup 
    3 except ImportError: 
    4     from distutils.core import setup 
     1import sys 
     2from setuptools import setup 
     3from setuptools.command.test import test as TestCommand 
     4 
     5class PyTest(TestCommand): 
     6    user_options = [('pytest-args=', 'a', "Arguments to pass to pytest")] 
     7 
     8    def initialize_options(self): 
     9        TestCommand.initialize_options(self) 
     10        self.pytest_args = '' 
     11 
     12    def run_tests(self): 
     13        import shlex 
     14        #import here, cause outside the eggs aren't loaded 
     15        import pytest 
     16        errno = pytest.main(shlex.split(self.pytest_args)) 
     17        sys.exit(errno) 
    518 
    619def find_version(package): 
    … …  
    4861        'sasmodels': ['*.c', '*.cl'], 
    4962    }, 
    50     install_requires = [ 
     63    install_requires=[ 
    5164    ], 
    52     extras_require = { 
     65    extras_require={ 
    5366        'OpenCL': ["pyopencl"], 
    5467        'Bumps': ["bumps"], 
    5568        'TinyCC': ["tinycc"], 
    5669    }, 
     70    build_requires=['setuptools'], 
     71    test_requires=['pytest'], 
     72    cmdclass = {'test': PyTest}, 
    5773) 

  • doc/guide/pd/polydispersity.rst

    reda8b30 r92d330fd  
    4242calculations are generally more robust with more data points or more angles. 
    4343 
    44 The following five distribution functions are provided: 
     44The following distribution functions are provided: 
    4545 
    4646*  *Rectangular Distribution* 
     47*  *Uniform Distribution* 
    4748*  *Gaussian Distribution* 
    4849*  *Lognormal Distribution* 
    4950*  *Schulz Distribution* 
    5051*  *Array Distribution* 
     52*  *Boltzmann Distribution* 
    5153 
    5254These are all implemented as *number-average* distributions. 
    … …  
    8587    Rectangular distribution. 
    8688 
     89 
     90 
     91Uniform Distribution 
     92^^^^^^^^^^^^^^^^^^^^^^^^ 
     93 
     94The Uniform Distribution is defined as 
     95 
     96    .. math:: 
     97 
     98        f(x) = \frac{1}{\text{Norm}} 
     99        \begin{cases} 
     100          1 & \text{for } |x - \bar x| \leq \sigma \\ 
     101          0 & \text{for } |x - \bar x| > \sigma 
     102        \end{cases} 
     103 
     104    where $\bar x$ is the mean of the distribution, $\sigma$ is the half-width, and 
     105    *Norm* is a normalization factor which is determined during the numerical 
     106    calculation. 
     107 
     108    Note that the polydispersity is given by 
     109 
     110    .. math:: \text{PD} = \sigma / \bar x 
     111 
     112    .. figure:: pd_uniform.jpg 
     113 
     114        Uniform distribution. 
     115 
     116The value $N_\sigma$ is ignored for this distribution. 
     117 
    87118.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
    88119 
    … …  
    183214^^^^^^^^^^^^^^^^^^ 
    184215 
    185 This user-definable distribution should be given as as a simple ASCII text 
     216This user-definable distribution should be given as a simple ASCII text 
    186217file where the array is defined by two columns of numbers: $x$ and $f(x)$. 
    187218The $f(x)$ will be normalized to 1 during the computation. 
    … …  
    202233given for the model will have no affect, and will be ignored when computing 
    203234the average.  This means that any parameter with an array distribution will 
    204 not be fittable. 
     235not be fitable. 
     236 
     237.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     238 
     239Boltzmann Distribution 
     240^^^^^^^^^^^^^^^^^^^^^^ 
     241 
     242The Boltzmann Distribution is defined as 
     243 
     244.. math:: 
     245 
     246    f(x) = \frac{1}{\text{Norm}} 
     247           \exp\left(-\frac{ | x - \bar x | }{\sigma}\right) 
     248 
     249where $\bar x$ is the mean of the distribution and *Norm* is a normalization 
     250factor which is determined during the numerical calculation. 
     251The width is defined as 
     252 
     253.. math:: \sigma=\frac{k T}{E} 
     254 
     255which is the inverse Boltzmann factor, 
     256where $k$ is the Boltzmann constant, $T$ the temperature in Kelvin and $E$ a 
     257characteristic energy per particle. 
     258 
     259.. figure:: pd_boltzmann.jpg 
     260 
     261    Boltzmann distribution. 
    205262 
    206263.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 

  • explore/realspace.py

    r8fb2a94 r032646d  
    33import time 
    44from copy import copy 
     5import os 
     6import argparse 
     7from collections import OrderedDict 
    58 
    69import numpy as np 
    710from numpy import pi, radians, sin, cos, sqrt 
    8 from numpy.random import poisson, uniform 
     11from numpy.random import poisson, uniform, randn, rand 
    912from numpy.polynomial.legendre import leggauss 
    1013from scipy.integrate import simps 
    1114from scipy.special import j1 as J1 
     15 
     16try: 
     17    import numba 
     18    USE_NUMBA = True 
     19except ImportError: 
     20    USE_NUMBA = False 
    1221 
    1322# Definition of rotation matrices comes from wikipedia: 
    … …  
    8291        raise NotImplementedError() 
    8392 
     93    def dims(self): 
     94        # type: () -> float, float, float 
     95        raise NotImplementedError() 
     96 
    8497    def rotate(self, theta, phi, psi): 
    8598        self.rotation = rotation(theta, phi, psi) * self.rotation 
    … …  
    116129                     for s2 in shapes] 
    117130        self.r_max = max(distances + [s.r_max for s in shapes]) 
    118  
    119     def volume(self): 
    120         return sum(shape.volume() for shape in self.shapes) 
     131        self.volume = sum(shape.volume for shape in self.shapes) 
    121132 
    122133    def sample(self, density): 
    … …  
    133144        self._scale = np.array([a/2, b/2, c/2])[None, :] 
    134145        self.r_max = sqrt(a**2 + b**2 + c**2) 
    135  
    136     def volume(self): 
    137         return self.a*self.b*self.c 
     146        self.dims = a, b, c 
     147        self.volume = a*b*c 
    138148 
    139149    def sample(self, density): 
    140         num_points = poisson(density*self.a*self.b*self.c) 
     150        num_points = poisson(density*self.volume) 
    141151        points = self._scale*uniform(-1, 1, size=(num_points, 3)) 
    142152        values = self.value.repeat(points.shape[0]) 
    … …  
    152162        self._scale = np.array([ra, rb, length/2])[None, :] 
    153163        self.r_max = sqrt(4*max(ra, rb)**2 + length**2) 
    154  
    155     def volume(self): 
    156         return pi*self.ra*self.rb*self.length 
     164        self.dims = 2*ra, 2*rb, length 
     165        self.volume = pi*ra*rb*length 
    157166 
    158167    def sample(self, density): 
    159         # density of the bounding box 
     168        # randomly sample from a box of side length 2*r, excluding anything 
     169        # not in the cylinder 
    160170        num_points = poisson(density*4*self.ra*self.rb*self.length) 
    161171        points = uniform(-1, 1, size=(num_points, 3)) 
    … …  
    174184        self._scale = np.array([ra, rb, rc])[None, :] 
    175185        self.r_max = 2*max(ra, rb, rc) 
    176  
    177     def volume(self): 
    178         return 4*pi/3 * self.ra * self.rb * self.rc 
     186        self.dims = 2*ra, 2*rb, 2*rc 
     187        self.volume = 4*pi/3 * ra * rb * rc 
    179188 
    180189    def sample(self, density): 
    … …  
    194203        self.rotate(*orientation) 
    195204        self.shift(*center) 
     205        helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2) 
     206        total_radius = self.helix_radius + self.tube_radius 
    196207        self.helix_radius, self.helix_pitch = helix_radius, helix_pitch 
    197208        self.tube_radius, self.tube_length = tube_radius, tube_length 
    198         helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2) 
    199         self.r_max = sqrt((2*helix_radius + 2*tube_radius)*2 
    200                           + (helix_length + 2*tube_radius)**2) 
    201  
    202     def volume(self): 
     209        self.r_max = sqrt(4*total_radius + (helix_length + 2*tube_radius)**2) 
     210        self.dims = 2*total_radius, 2*total_radius, helix_length 
    203211        # small tube radius approximation; for larger tubes need to account 
    204212        # for the fact that the inner length is much shorter than the outer 
    205213        # length 
    206         return pi*self.tube_radius**2*self.tube_length 
     214        self.volume = pi*self.tube_radius**2*self.tube_length 
    207215 
    208216    def points(self, density): 
    … …  
    244252        return values, self._adjust(points) 
    245253 
    246 NUMBA = False 
    247 if NUMBA: 
     254def csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 
     255    core = Box(a, b, c, sld_core) 
     256    side_a = Box(da, b, c, slda, center=((a+da)/2, 0, 0)) 
     257    side_b = Box(a, db, c, sldb, center=(0, (b+db)/2, 0)) 
     258    side_c = Box(a, b, dc, sldc, center=(0, 0, (c+dc)/2)) 
     259    side_a2 = copy(side_a).shift(-a-da, 0, 0) 
     260    side_b2 = copy(side_b).shift(0, -b-db, 0) 
     261    side_c2 = copy(side_c).shift(0, 0, -c-dc) 
     262    shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 
     263    shape.dims = 2*da+a, 2*db+b, 2*dc+c 
     264    return shape 
     265 
     266def _Iqxy(values, x, y, z, qa, qb, qc): 
     267    """I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r)""" 
     268    Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2 
     269            for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)] 
     270    return Iq 
     271 
     272if USE_NUMBA: 
     273    # Override simple numpy solution with numba if available 
    248274    from numba import njit 
    249275    @njit("f8[:](f8[:],f8[:],f8[:],f8[:],f8[:],f8[:],f8[:])") 
    … …  
    252278        for j in range(len(Iq)): 
    253279            total = 0. + 0j 
    254             for k in range(len(Iq)): 
     280            for k in range(len(values)): 
    255281                total += values[k]*np.exp(1j*(qa[j]*x[k] + qb[j]*y[k] + qc[j]*z[k])) 
    256282            Iq[j] = abs(total)**2 
    … …  
    263289    values = rho*volume 
    264290    x, y, z = points.T 
     291    values, x, y, z, qa, qb, qc = [np.asarray(v, 'd') 
     292                                   for v in (values, x, y, z, qa, qb, qc)] 
    265293 
    266294    # I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r) 
    267     if NUMBA: 
    268         Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten()) 
    269     else: 
    270         Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2 
    271               for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)] 
     295    Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten()) 
    272296    return np.asarray(Iq).reshape(qx.shape) / np.sum(volume) 
    273297 
    … …  
    338362""" 
    339363 
    340 if NUMBA: 
     364if USE_NUMBA: 
     365    # Override simple numpy solution with numba if available 
    341366    @njit("f8[:](f8[:], f8[:], f8[:,:])") 
    342     def _calc_Pr_uniform_jit(r, rho, points): 
     367    def _calc_Pr_uniform(r, rho, points): 
    343368        dr = r[0] 
    344369        n_max = len(r) 
    … …  
    362387    # P(r) with uniform steps in r is 3x faster; check if we are uniform 
    363388    # before continuing 
     389    r, rho, points = [np.asarray(v, 'd') for v in (r, rho, points)] 
    364390    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01: 
    365391        Pr = _calc_Pr_nonuniform(r, rho, points) 
    366392    else: 
    367         if NUMBA: 
    368             Pr = _calc_Pr_uniform_jit(r, rho, points) 
    369         else: 
    370             Pr = _calc_Pr_uniform(r, rho, points) 
     393        Pr = _calc_Pr_uniform(r, rho, points) 
    371394    return Pr / Pr.max() 
    372395 
    … …  
    399422    return edges 
    400423 
     424# -------------- plotters ---------------- 
    401425def plot_calc(r, Pr, q, Iq, theory=None): 
    402426    import matplotlib.pyplot as plt 
    … …  
    410434    plt.ylabel('Iq') 
    411435    if theory is not None: 
    412         plt.loglog(theory[0], theory[1], '-', label='analytic') 
     436        plt.loglog(theory[0], theory[1]/theory[1][0], '-', label='analytic') 
    413437        plt.legend() 
    414438 
    … …  
    444468    ax.autoscale(True) 
    445469 
    446 def check_shape(shape, fn=None): 
    447     rho_solvent = 0 
    448     q = np.logspace(-3, 0, 200) 
    449     r = shape.r_bins(q, r_step=0.01) 
    450     sampling_density = 6*5000 / shape.volume() 
    451     rho, points = shape.sample(sampling_density) 
    452     t0 = time.time() 
    453     Pr = calc_Pr(r, rho-rho_solvent, points) 
    454     print("calc Pr time", time.time() - t0) 
    455     Iq = calc_Iq(q, r, Pr) 
    456     theory = (q, fn(q)) if fn is not None else None 
    457  
    458     import pylab 
    459     #plot_points(rho, points); pylab.figure() 
    460     plot_calc(r, Pr, q, Iq, theory=theory) 
    461     pylab.show() 
    462  
    463 def check_shape_2d(shape, fn=None, view=(0, 0, 0)): 
    464     rho_solvent = 0 
    465     nq, qmax = 100, 1.0 
    466     qx = np.linspace(0.0, qmax, nq) 
    467     qy = np.linspace(0.0, qmax, nq) 
    468     Qx, Qy = np.meshgrid(qx, qy) 
    469     sampling_density = 50000 / shape.volume() 
    470     #t0 = time.time() 
    471     rho, points = shape.sample(sampling_density) 
    472     #print("sample time", time.time() - t0) 
    473     t0 = time.time() 
    474     Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 
    475     print("calc time", time.time() - t0) 
    476     theory = fn(Qx, Qy) if fn is not None else None 
    477     Iqxy += 0.001 * Iqxy.max() 
    478     if theory is not None: 
    479         theory += 0.001 * theory.max() 
    480  
    481     import pylab 
    482     #plot_points(rho, points); pylab.figure() 
    483     plot_calc_2d(qx, qy, Iqxy, theory=theory) 
    484     pylab.show() 
    485  
     470# ----------- Analytic models -------------- 
    486471def sas_sinx_x(x): 
    487472    with np.errstate(all='ignore'): 
    … …  
    510495    for k, qk in enumerate(q): 
    511496        qab, qc = qk*sin_alpha, qk*cos_alpha 
    512         Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2) 
     497        Fq = sas_2J1x_x(qab*radius) * sas_sinx_x(qc*length/2) 
    513498        Iq[k] = np.sum(w*Fq**2) 
    514     Iq = Iq/Iq[0] 
     499    Iq = Iq 
    515500    return Iq 
    516501 
    517502def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)): 
    518503    qa, qb, qc = invert_view(qx, qy, view) 
    519     qab = np.sqrt(qa**2 + qb**2) 
    520     Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2) 
     504    qab = sqrt(qa**2 + qb**2) 
     505    Fq = sas_2J1x_x(qab*radius) * sas_sinx_x(qc*length/2) 
    521506    Iq = Fq**2 
    522507    return Iq.reshape(qx.shape) 
    … …  
    524509def sphere_Iq(q, radius): 
    525510    Iq = sas_3j1x_x(q*radius)**2 
    526     return Iq/Iq[0] 
     511    return Iq 
     512 
     513def box_Iq(q, a, b, c): 
     514    z, w = leggauss(76) 
     515    outer_sum = np.zeros_like(q) 
     516    for cos_alpha, outer_w in zip((z+1)/2, w): 
     517        sin_alpha = sqrt(1.0-cos_alpha*cos_alpha) 
     518        qc = q*cos_alpha 
     519        siC = c*sas_sinx_x(c*qc/2) 
     520        inner_sum = np.zeros_like(q) 
     521        for beta, inner_w in zip((z + 1)*pi/4, w): 
     522            qa, qb = q*sin_alpha*sin(beta), q*sin_alpha*cos(beta) 
     523            siA = a*sas_sinx_x(a*qa/2) 
     524            siB = b*sas_sinx_x(b*qb/2) 
     525            Fq = siA*siB*siC 
     526            inner_sum += inner_w * Fq**2 
     527        outer_sum += outer_w * inner_sum 
     528    Iq = outer_sum / 4  # = outer*um*zm*8.0/(4.0*M_PI) 
     529    return Iq 
     530 
     531def box_Iqxy(qx, qy, a, b, c, view=(0, 0, 0)): 
     532    qa, qb, qc = invert_view(qx, qy, view) 
     533    sia = sas_sinx_x(qa*a/2) 
     534    sib = sas_sinx_x(qb*b/2) 
     535    sic = sas_sinx_x(qc*c/2) 
     536    Fq = sia*sib*sic 
     537    Iq = Fq**2 
     538    return Iq.reshape(qx.shape) 
    527539 
    528540def csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core): 
    … …  
    539551        sin_alpha = sqrt(1.0-cos_alpha*cos_alpha) 
    540552        qc = q*cos_alpha 
    541         siC = c*j0(c*qc/2) 
    542         siCt = tC*j0(tC*qc/2) 
     553        siC = c*sas_sinx_x(c*qc/2) 
     554        siCt = tC*sas_sinx_x(tC*qc/2) 
    543555        inner_sum = np.zeros_like(q) 
    544556        for beta, inner_w in zip((z + 1)*pi/4, w): 
    545557            qa, qb = q*sin_alpha*sin(beta), q*sin_alpha*cos(beta) 
    546             siA = a*j0(a*qa/2) 
    547             siB = b*j0(b*qb/2) 
    548             siAt = tA*j0(tA*qa/2) 
    549             siBt = tB*j0(tB*qb/2) 
     558            siA = a*sas_sinx_x(a*qa/2) 
     559            siB = b*sas_sinx_x(b*qb/2) 
     560            siAt = tA*sas_sinx_x(tA*qa/2) 
     561            siBt = tB*sas_sinx_x(tB*qb/2) 
    550562            if overlapping: 
    551563                Fq = (dr0*siA*siB*siC 
    … …  
    584596    return Iq.reshape(qx.shape) 
    585597 
    586 def check_cylinder(radius=25, length=125, rho=2.): 
     598# --------- Test cases ----------- 
     599 
     600def build_cylinder(radius=25, length=125, rho=2.): 
    587601    shape = EllipticalCylinder(radius, radius, length, rho) 
    588     fn = lambda q: cylinder_Iq(q, radius, length) 
    589     check_shape(shape, fn) 
    590  
    591 def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)): 
    592     shape = EllipticalCylinder(radius, radius, length, rho) 
    593     fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
    594     check_shape_2d(shape, fn, view=view) 
    595  
    596 def check_cylinder_2d_lattice(radius=25, length=125, rho=2., 
    597                               view=(0, 0, 0)): 
    598     nx, dx = 1, 2*radius 
    599     ny, dy = 30, 2*radius 
    600     nz, dz = 30, length 
    601     dx, dy, dz = 2*dx, 2*dy, 2*dz 
    602     def center(*args): 
    603         sigma = 0.333 
    604         space = 2 
    605         return [(space*n+np.random.randn()*sigma)*x for n, x in args] 
    606     shapes = [EllipticalCylinder(radius, radius, length, rho, 
    607                                  #center=(ix*dx, iy*dy, iz*dz) 
    608                                  orientation=np.random.randn(3)*0, 
    609                                  center=center((ix, dx), (iy, dy), (iz, dz)) 
    610                                 ) 
     602    fn = lambda q: cylinder_Iq(q, radius, length)*rho**2 
     603    fn_xy = lambda qx, qy, view: cylinder_Iqxy(qx, qy, radius, length, view=view)*rho**2 
     604    return shape, fn, fn_xy 
     605 
     606def build_sphere(radius=125, rho=2): 
     607    shape = TriaxialEllipsoid(radius, radius, radius, rho) 
     608    fn = lambda q: sphere_Iq(q, radius)*rho**2 
     609    fn_xy = lambda qx, qy, view: sphere_Iq(np.sqrt(qx**2+qy**2), radius)*rho**2 
     610    return shape, fn, fn_xy 
     611 
     612def build_box(a=10, b=20, c=30, rho=2.): 
     613    shape = Box(a, b, c, rho) 
     614    fn = lambda q: box_Iq(q, a, b, c)*rho**2 
     615    fn_xy = lambda qx, qy, view: box_Iqxy(qx, qy, a, b, c, view=view)*rho**2 
     616    return shape, fn, fn_xy 
     617 
     618def build_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 
     619    shape = csbox(a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
     620    fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
     621    fn_xy = lambda qx, qy, view: csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 
     622                                            slda, sldb, sldc, sld_core, view=view) 
     623    return shape, fn, fn_xy 
     624 
     625def build_cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     626                  shuffle=0, rotate=0): 
     627    a, b, c = shape.dims 
     628    shapes = [copy(shape) 
     629              .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     630                     (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     631                     (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     632              .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
    611633              for ix in range(nx) 
    612634              for iy in range(ny) 
    613635              for iz in range(nz)] 
    614     shape = Composite(shapes) 
    615     fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
    616     check_shape_2d(shape, fn, view=view) 
    617  
    618 def check_sphere(radius=125, rho=2): 
    619     shape = TriaxialEllipsoid(radius, radius, radius, rho) 
    620     fn = lambda q: sphere_Iq(q, radius) 
    621     check_shape(shape, fn) 
    622  
    623 def check_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 
    624     core = Box(a, b, c, sld_core) 
    625     side_a = Box(da, b, c, slda, center=((a+da)/2, 0, 0)) 
    626     side_b = Box(a, db, c, sldb, center=(0, (b+db)/2, 0)) 
    627     side_c = Box(a, b, dc, sldc, center=(0, 0, (c+dc)/2)) 
    628     side_a2 = copy(side_a).shift(-a-da, 0, 0) 
    629     side_b2 = copy(side_b).shift(0, -b-db, 0) 
    630     side_c2 = copy(side_c).shift(0, 0, -c-dc) 
    631     shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 
    632     def fn(q): 
    633         return csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
    634     #check_shape(shape, fn) 
    635  
    636     view = (20, 30, 40) 
    637     def fn_xy(qx, qy): 
    638         return csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 
    639                           slda, sldb, sldc, sld_core, view=view) 
    640     check_shape_2d(shape, fn_xy, view=view) 
     636    lattice = Composite(shapes) 
     637    return lattice 
     638 
     639 
     640SHAPE_FUNCTIONS = OrderedDict([ 
     641    ("cylinder", build_cylinder), 
     642    ("sphere", build_sphere), 
     643    ("box", build_box), 
     644    ("csbox", build_csbox), 
     645]) 
     646SHAPES = list(SHAPE_FUNCTIONS.keys()) 
     647 
     648def check_shape(title, shape, fn=None, show_points=False, 
     649                mesh=100, qmax=1.0, r_step=0.01, samples=5000): 
     650    rho_solvent = 0 
     651    qmin = qmax/100. 
     652    q = np.logspace(np.log10(qmin), np.log10(qmax), mesh) 
     653    r = shape.r_bins(q, r_step=r_step) 
     654    sampling_density = samples / shape.volume 
     655    rho, points = shape.sample(sampling_density) 
     656    t0 = time.time() 
     657    Pr = calc_Pr(r, rho-rho_solvent, points) 
     658    print("calc Pr time", time.time() - t0) 
     659    Iq = calc_Iq(q, r, Pr) 
     660    theory = (q, fn(q)) if fn is not None else None 
     661 
     662    import pylab 
     663    if show_points: 
     664         plot_points(rho, points); pylab.figure() 
     665    plot_calc(r, Pr, q, Iq, theory=theory) 
     666    pylab.gcf().canvas.set_window_title(title) 
     667    pylab.show() 
     668 
     669def check_shape_2d(title, shape, fn=None, view=(0, 0, 0), show_points=False, 
     670                   mesh=100, qmax=1.0, samples=5000): 
     671    rho_solvent = 0 
     672    qx = np.linspace(0.0, qmax, mesh) 
     673    qy = np.linspace(0.0, qmax, mesh) 
     674    Qx, Qy = np.meshgrid(qx, qy) 
     675    sampling_density = samples / shape.volume 
     676    t0 = time.time() 
     677    rho, points = shape.sample(sampling_density) 
     678    print("point generation time", time.time() - t0) 
     679    t0 = time.time() 
     680    Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 
     681    print("calc Iqxy time", time.time() - t0) 
     682    theory = fn(Qx, Qy, view) if fn is not None else None 
     683    Iqxy += 0.001 * Iqxy.max() 
     684    if theory is not None: 
     685        theory += 0.001 * theory.max() 
     686 
     687    import pylab 
     688    if show_points: 
     689        plot_points(rho, points); pylab.figure() 
     690    plot_calc_2d(qx, qy, Iqxy, theory=theory) 
     691    pylab.gcf().canvas.set_window_title(title) 
     692    pylab.show() 
     693 
     694def main(): 
     695    parser = argparse.ArgumentParser( 
     696        description="Compute scattering from realspace sampling", 
     697        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
     698        ) 
     699    parser.add_argument('-d', '--dim', type=int, default=1, 
     700                        help='dimension 1 or 2') 
     701    parser.add_argument('-m', '--mesh', type=int, default=100, 
     702                        help='number of mesh points') 
     703    parser.add_argument('-s', '--samples', type=int, default=5000, 
     704                        help="number of sample points") 
     705    parser.add_argument('-q', '--qmax', type=float, default=0.5, 
     706                        help='max q') 
     707    parser.add_argument('-v', '--view', type=str, default='0,0,0', 
     708                        help='theta,phi,psi angles') 
     709    parser.add_argument('-n', '--lattice', type=str, default='1,1,1', 
     710                        help='lattice size') 
     711    parser.add_argument('-z', '--spacing', type=str, default='2,2,2', 
     712                        help='lattice spacing') 
     713    parser.add_argument('-r', '--rotate', type=float, default=0., 
     714                        help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise") 
     715    parser.add_argument('-w', '--shuffle', type=float, default=0., 
     716                        help="position relative to lattice, gaussian < 0.3, uniform otherwise") 
     717    parser.add_argument('-p', '--plot', action='store_true', 
     718                        help='plot points') 
     719    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], 
     720                        help='oriented shape') 
     721    parser.add_argument('pars', type=str, nargs='*', help='shape parameters') 
     722    opts = parser.parse_args() 
     723    pars = {key: float(value) for p in opts.pars for key, value in [p.split('=')]} 
     724    nx, ny, nz = [int(v) for v in opts.lattice.split(',')] 
     725    dx, dy, dz = [float(v) for v in opts.spacing.split(',')] 
     726    shuffle, rotate = opts.shuffle, opts.rotate 
     727    shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 
     728    if nx > 1 or ny > 1 or nz > 1: 
     729        shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 
     730    title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 
     731    if opts.dim == 1: 
     732        check_shape(title, shape, fn, show_points=opts.plot, 
     733                    mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 
     734    else: 
     735        view = tuple(float(v) for v in opts.view.split(',')) 
     736        check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot, 
     737                       mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 
     738 
    641739 
    642740if __name__ == "__main__": 
    643     check_cylinder(radius=10, length=40) 
    644     #check_cylinder_2d(radius=10, length=40, view=(90,30,0)) 
    645     #check_cylinder_2d_lattice(radius=10, length=50, view=(90,30,0)) 
    646     #check_sphere() 
    647     #check_csbox() 
    648     #check_csbox(da=100, db=200, dc=300) 
     741    main() 

  • sasmodels/kernel_iq.c

    r924a119 raadec17  
    9393// Compute the magnetic sld 
    9494static double mag_sld( 
    95   const int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 
     95  const unsigned int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 
    9696  const double qx, const double qy, 
    9797  const double px, const double py, 
    … …  
    103103    const double perp = qy*mx - qx*my; 
    104104    switch (xs) { 
     105      default: // keep compiler happy; condition ensures xs in [0,1,2,3] 
    105106      case 0: // uu => sld - D M_perpx 
    106107          return sld - px*perp; 
    … …  
    659660 
    660661            // loop over uu, ud real, du real, dd, ud imag, du imag 
    661             for (int xs=0; xs<6; xs++) { 
     662            for (unsigned int xs=0; xs<6; xs++) { 
    662663              const double xs_weight = spins[xs]; 
    663664              if (xs_weight > 1.e-8) { 

  • sasmodels/weights.py

    r2d81cfe r3d58247  
    9797        return x, px 
    9898 
     99class UniformDispersion(Dispersion): 
     100    r""" 
     101    Uniform dispersion, with width $\sigma$. 
     102 
     103    .. math:: 
     104 
     105        w = 1 
     106    """ 
     107    type = "uniform" 
     108    default = dict(npts=35, width=0, nsigmas=None) 
     109    def _weights(self, center, sigma, lb, ub): 
     110        x = np.linspace(center-sigma, center+sigma, self.npts) 
     111        x = x[(x >= lb) & (x <= ub)] 
     112        return x, np.ones_like(x) 
    99113 
    100114class RectangleDispersion(Dispersion): 
    … …  
    107121    """ 
    108122    type = "rectangle" 
    109     default = dict(npts=35, width=0, nsigmas=1.70325) 
    110     def _weights(self, center, sigma, lb, ub): 
    111         x = self._linspace(center, sigma, lb, ub) 
    112         x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)] 
    113         return x, np.ones_like(x) 
    114  
     123    default = dict(npts=35, width=0, nsigmas=1.73205) 
     124    def _weights(self, center, sigma, lb, ub): 
     125         x = self._linspace(center, sigma, lb, ub) 
     126         x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)] 
     127         return x, np.ones_like(x) 
    115128 
    116129class LogNormalDispersion(Dispersion): 
    … …  
    190203        return x, px 
    191204 
     205class BoltzmannDispersion(Dispersion): 
     206    r""" 
     207    Boltzmann dispersion, with $\sigma=k T/E$. 
     208 
     209    .. math:: 
     210 
     211        w = \exp\left( -|x - c|/\sigma\right) 
     212    """ 
     213    type = "boltzmann" 
     214    default = dict(npts=35, width=0, nsigmas=3) 
     215    def _weights(self, center, sigma, lb, ub): 
     216        x = self._linspace(center, sigma, lb, ub) 
     217        px = np.exp(-np.abs(x-center) / np.abs(sigma)) 
     218        return x, px 
    192219 
    193220# dispersion name -> disperser lookup table. 
    … …  
    196223MODELS = OrderedDict((d.type, d) for d in ( 
    197224    RectangleDispersion, 
     225    UniformDispersion, 
    198226    ArrayDispersion, 
    199227    LogNormalDispersion, 
    200228    GaussianDispersion, 
    201229    SchulzDispersion, 
     230    BoltzmannDispersion 
    202231)) 
    203232 
Note: See TracChangeset for help on using the changeset viewer.