Changeset 108e70e in sasmodels for sasmodels


Ignore:
Timestamp:
Dec 14, 2017 1:08:45 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
ef85a09
Parents:
df69efa
Message:

use Iqac/Iqabc? for the new orientation interface but Iqxy for the old

Location:
sasmodels
Files:
29 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/details.py

    r2d81cfe r108e70e  
    258258    # type: (...) -> Sequence[np.ndarray] 
    259259    """ 
     260    **Deprecated** Theta weights will be computed in the kernel wrapper if 
     261    they are needed. 
     262 
    260263    If there is a theta parameter, update the weights of that parameter so that 
    261264    the cosine weighting required for polar integration is preserved. 
     
    272275    Returns updated weights vectors 
    273276    """ 
    274     # TODO: explain in a comment why scale and background are missing 
    275277    # Apparently the parameters.theta_offset similarly skips scale and 
    276278    # and background, so the indexing works out, but they are still shipped 
     
    279281        index = parameters.theta_offset 
    280282        theta = dispersity[index] 
    281         # TODO: modify the dispersity vector to avoid the theta=-90,90,270,... 
    282283        theta_weight = abs(cos(radians(theta))) 
    283284        weights = tuple(theta_weight*w if k == index else w 
  • sasmodels/generate.py

    ra261a83 r108e70e  
    77    particular dimensions averaged over all orientations. 
    88 
    9     *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form 
    10     with particular dimensions for a single orientation. 
    11  
    12     *Imagnetic(qx, qy, result[], p1, p2, ...)* returns the scattering for the 
    13     polarized neutron spin states (up-up, up-down, down-up, down-down) for 
    14     a form with particular dimensions for a single orientation. 
     9    *Iqac(qab, qc, p1, p2, ...)* returns the scattering at qab, qc 
     10    for a rotationally symmetric form with particular dimensions. 
     11    qab, qc are determined from shape orientation and scattering angles. 
     12    This call is used if the shape has orientation parameters theta and phi. 
     13 
     14    *Iqabc(qa, qb, qc, p1, p2, ...)* returns the scattering at qa, qb, qc 
     15    for a form with particular dimensions.  qa, qb, qc are determined from 
     16    shape orientation and scattering angles. This call is used if the shape 
     17    has orientation parameters theta, phi and psi. 
     18 
     19    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy.  Use this 
     20    to create an arbitrary 2D theory function, needed for q-dependent 
     21    background functions and for models with non-uniform magnetism. 
    1522 
    1623    *form_volume(p1, p2, ...)* returns the volume of the form with particular 
     
    3138scale and background parameters for each model. 
    3239 
    33 *Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99 
    34 functions written for OpenCL.  All functions need prototype declarations 
    35 even if the are defined before they are used.  OpenCL does not support 
    36 *#include* preprocessor directives, so instead the list of includes needs 
    37 to be given as part of the metadata in the kernel module definition. 
    38 The included files should be listed using a path relative to the kernel 
    39 module, or if using "lib/file.c" if it is one of the standard includes 
    40 provided with the sasmodels source.  The includes need to be listed in 
    41 order so that functions are defined before they are used. 
     40C code should be stylized C-99 functions written for OpenCL. All functions 
     41need prototype declarations even if the are defined before they are used. 
     42Although OpenCL supports *#include* preprocessor directives, the list of 
     43includes should be given as part of the metadata in the kernel module 
     44definition. The included files should be listed using a path relative to the 
     45kernel module, or if using "lib/file.c" if it is one of the standard includes 
     46provided with the sasmodels source. The includes need to be listed in order 
     47so that functions are defined before they are used. 
    4248 
    4349Floating point values should be declared as *double*.  For single precision 
     
    107113    present, the volume ratio is 1. 
    108114 
    109     *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the 
    110     C source code for the body of the volume, Iq, and Iqxy functions 
     115    *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing 
     116    the C source code for the body of the volume, Iq, and Iqac functions 
    111117    respectively.  These can also be defined in the last source file. 
    112118 
    113     *Iq* and *Iqxy* also be instead be python functions defining the 
     119    *Iq*, *Iqac*, *Iqabc* also be instead be python functions defining the 
    114120    kernel.  If they are marked as *Iq.vectorized = True* then the 
    115121    kernel is passed the entire *q* vector at once, otherwise it is 
     
    168174from zlib import crc32 
    169175from inspect import currentframe, getframeinfo 
     176import logging 
    170177 
    171178import numpy as np  # type: ignore 
     
    181188    pass 
    182189# pylint: enable=unused-import 
     190 
     191logger = logging.getLogger(__name__) 
    183192 
    184193# jitter projection to use in the kernel code.  See explore/jitter.py 
     
    626635 
    627636""" 
    628 def _gen_fn(name, pars, body, filename, line): 
    629     # type: (str, List[Parameter], str, str, int) -> str 
     637def _gen_fn(model_info, name, pars): 
     638    # type: (ModelInfo, str, List[Parameter]) -> str 
    630639    """ 
    631640    Generate a function given pars and body. 
     
    639648    """ 
    640649    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 
     650    body = getattr(model_info, name) 
     651    filename = model_info.filename 
     652    # Note: if symbol is defined strangely in the module then default it to 1 
     653    lineno = model_info.lineno.get(name, 1) 
    641654    return _FN_TEMPLATE % { 
    642655        'name': name, 'pars': par_decl, 'body': body, 
    643         'filename': filename.replace('\\', '\\\\'), 'line': line, 
     656        'filename': filename.replace('\\', '\\\\'), 'line': lineno, 
    644657    } 
    645658 
     
    656669 
    657670# type in IQXY pattern could be single, float, double, long double, ... 
    658 _IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)", 
     671_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
    659672                           flags=re.MULTILINE) 
    660 def _have_Iqxy(sources): 
     673def find_xy_mode(source): 
    661674    # type: (List[str]) -> bool 
    662675    """ 
    663     Return true if any file defines Iqxy. 
     676    Return the xy mode as qa, qac, qabc or qxy. 
    664677 
    665678    Note this is not a C parser, and so can be easily confused by 
    666679    non-standard syntax.  Also, it will incorrectly identify the following 
    667     as having Iqxy:: 
     680    as having 2D models:: 
    668681 
    669682        /* 
    670         double Iqxy(qx, qy, ...) { ... fill this in later ... } 
     683        double Iqac(qab, qc, ...) { ... fill this in later ... } 
    671684        */ 
    672685 
    673     If you want to comment out an Iqxy function, use // on the front of the 
    674     line instead. 
    675     """ 
    676     for _path, code in sources: 
    677         if _IQXY_PATTERN.search(code): 
    678             return True 
    679     return False 
    680  
    681  
    682 def _add_source(source, code, path): 
     686    If you want to comment out the function, use // on the front of the 
     687    line:: 
     688 
     689        /* 
     690        // double Iqac(qab, qc, ...) { ... fill this in later ... } 
     691        */ 
     692 
     693    """ 
     694    for code in source: 
     695        m = _IQXY_PATTERN.search(code) 
     696        if m is not None: 
     697            return m.group('mode') 
     698    return 'qa' 
     699 
     700 
     701def _add_source(source, code, path, lineno=1): 
    683702    """ 
    684703    Add a file to the list of source code chunks, tagged with path and line. 
    685704    """ 
    686705    path = path.replace('\\', '\\\\') 
    687     source.append('#line 1 "%s"' % path) 
     706    source.append('#line %d "%s"' % (lineno, path)) 
    688707    source.append(code) 
    689708 
     
    716735    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
    717736 
    718     # What kind of 2D model do we need? 
    719     xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str) 
    720                else 'qac' if not partable.is_asymmetric 
    721                else 'qabc') 
    722  
    723737    # Build initial sources 
    724738    source = [] 
     
    727741        _add_source(source, code, path) 
    728742 
     743    if model_info.c_code: 
     744        _add_source(source, model_info.c_code, model_info.filename, 
     745                    lineno=model_info.lineno.get('c_code', 1)) 
     746 
    729747    # Make parameters for q, qx, qy so that we can use them in declarations 
    730     q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
     748    q, qx, qy, qab, qa, qb, qc \ 
     749        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
    731750    # Generate form_volume function, etc. from body only 
    732751    if isinstance(model_info.form_volume, str): 
    733752        pars = partable.form_volume_parameters 
    734         source.append(_gen_fn('form_volume', pars, model_info.form_volume, 
    735                               model_info.filename, model_info._form_volume_line)) 
     753        source.append(_gen_fn(model_info, 'form_volume', pars)) 
    736754    if isinstance(model_info.Iq, str): 
    737755        pars = [q] + partable.iq_parameters 
    738         source.append(_gen_fn('Iq', pars, model_info.Iq, 
    739                               model_info.filename, model_info._Iq_line)) 
     756        source.append(_gen_fn(model_info, 'Iq', pars)) 
    740757    if isinstance(model_info.Iqxy, str): 
    741         pars = [qx, qy] + partable.iqxy_parameters 
    742         source.append(_gen_fn('Iqxy', pars, model_info.Iqxy, 
    743                               model_info.filename, model_info._Iqxy_line)) 
     758        pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters 
     759        source.append(_gen_fn(model_info, 'Iqxy', pars)) 
     760    if isinstance(model_info.Iqac, str): 
     761        pars = [qab, qc] + partable.iq_parameters 
     762        source.append(_gen_fn(model_info, 'Iqac', pars)) 
     763    if isinstance(model_info.Iqabc, str): 
     764        pars = [qa, qb, qc] + partable.iq_parameters 
     765        source.append(_gen_fn(model_info, 'Iqabc', pars)) 
     766 
     767    # What kind of 2D model do we need?  Is it consistent with the parameters? 
     768    xy_mode = find_xy_mode(source) 
     769    if xy_mode == 'qabc' and not partable.is_asymmetric: 
     770        raise ValueError("asymmetric oriented models need to define Iqabc") 
     771    elif xy_mode == 'qac' and partable.is_asymmetric: 
     772        raise ValueError("symmetric oriented models need to define Iqac") 
     773    elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'): 
     774        raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode) 
     775    elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'): 
     776        if xy_mode == 'qxy': 
     777            logger.warn("oriented shapes should define Iqac or Iqabc") 
     778        else: 
     779            raise ValueError("Expected function Iqac or Iqabc for oriented shape") 
    744780 
    745781    # Define the parameter table 
     
    767803    if xy_mode == 'qabc': 
    768804        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    769         call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
     805        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars 
    770806        clear_iqxy = "#undef CALL_IQ_ABC" 
    771807    elif xy_mode == 'qac': 
    772808        pars = ",".join(["_qa", "_qc"] + model_refs) 
    773         call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 
     809        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    774810        clear_iqxy = "#undef CALL_IQ_AC" 
    775     else:  # xy_mode == 'qa' 
     811    elif xy_mode == 'qa': 
    776812        pars = ",".join(["_qa"] + model_refs) 
    777813        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    778814        clear_iqxy = "#undef CALL_IQ_A" 
     815    elif xy_mode == 'qxy': 
     816        orientation_refs = _call_pars("_v.", partable.orientation_parameters) 
     817        pars = ",".join(["_qx", "_qy"] + model_refs + orientation_refs) 
     818        call_iqxy = "#define CALL_IQ_XY(_qx,_qy,_v) Iqxy(%s)" % pars 
     819        clear_iqxy = "#undef CALL_IQ_XY" 
     820        if partable.orientation_parameters: 
     821            call_iqxy += "\n#define HAVE_THETA" 
     822            clear_iqxy += "\n#undef HAVE_THETA" 
     823        if partable.is_asymmetric: 
     824            call_iqxy += "\n#define HAVE_PSI" 
     825            clear_iqxy += "\n#undef HAVE_PSI" 
     826 
    779827 
    780828    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
  • sasmodels/kernel_header.c

    r8698a0d r108e70e  
    150150inline double cube(double x) { return x*x*x; } 
    151151inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 
     152 
     153// CRUFT: support old style models with orientation received qx, qy and angles 
     154 
     155// To rotate from the canonical position to theta, phi, psi, first rotate by 
     156// psi about the major axis, oriented along z, which is a rotation in the 
     157// detector plane xy. Next rotate by theta about the y axis, aligning the major 
     158// axis in the xz plane. Finally, rotate by phi in the detector plane xy. 
     159// To compute the scattering, undo these rotations in reverse order: 
     160//     rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi 
     161// The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit 
     162// vector in the q direction. 
     163// To change between counterclockwise and clockwise rotation, change the 
     164// sign of phi and psi. 
     165 
     166#if 1 
     167//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
     168#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
     169    SINCOS(phi*M_PI_180, sn, cn); \ 
     170    q = sqrt(qx*qx + qy*qy); \ 
     171    cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
     172    sn = sqrt(1 - cn*cn); \ 
     173    } while (0) 
     174#else 
     175// SasView 3.x definition of orientation 
     176#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
     177    SINCOS(theta*M_PI_180, sn, cn); \ 
     178    q = sqrt(qx*qx + qy*qy);\ 
     179    cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \ 
     180    sn = sqrt(1 - cn*cn); \ 
     181    } while (0) 
     182#endif 
     183 
     184#if 1 
     185#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \ 
     186    q = sqrt(qx*qx + qy*qy); \ 
     187    const double qxhat = qx/q; \ 
     188    const double qyhat = qy/q; \ 
     189    double sin_theta, cos_theta; \ 
     190    double sin_phi, cos_phi; \ 
     191    double sin_psi, cos_psi; \ 
     192    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
     193    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
     194    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
     195    xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \ 
     196         + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \ 
     197    yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \ 
     198         + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \ 
     199    zhat = qxhat*(-sin_theta*cos_phi) \ 
     200         + qyhat*(-sin_theta*sin_phi); \ 
     201    } while (0) 
     202#else 
     203// SasView 3.x definition of orientation 
     204#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
     205    q = sqrt(qx*qx + qy*qy); \ 
     206    const double qxhat = qx/q; \ 
     207    const double qyhat = qy/q; \ 
     208    double sin_theta, cos_theta; \ 
     209    double sin_phi, cos_phi; \ 
     210    double sin_psi, cos_psi; \ 
     211    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
     212    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
     213    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
     214    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ 
     215    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ 
     216    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
     217    } while (0) 
     218#endif 
  • sasmodels/kernel_iq.c

    r6aee3ab r108e70e  
    3131//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 
    3232//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes 
     33//  CALL_IQ_XY(qx, qy, table) : call the Iqxy function for arbitrary models 
    3334//  INVALID(table) : test if the current point is feesible to calculate.  This 
    3435//      will be defined in the kernel definition file. 
     
    469470  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
    470471  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     472#elif defined(CALL_IQ_XY) 
     473  // direct call to qx,qy calculator 
     474  double qx, qy; 
     475  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     476  #define BUILD_ROTATION() do {} while(0) 
     477  #define APPLY_ROTATION() do {} while(0) 
     478  #define CALL_KERNEL() CALL_IQ_XY(qx, qy, local_values.table) 
    471479#endif 
    472480 
     
    477485#if defined(CALL_IQ) || defined(CALL_IQ_A) 
    478486  #define APPLY_PROJECTION() const double weight=weight0 
     487#elif defined(CALL_IQ_XY) 
     488  // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc 
     489  // Need to plug the values for the orientation angles back into parameter 
     490  // table in case they were overridden by the orientation offset.  This 
     491  // means that orientation dispersity will not work for these models, but 
     492  // it was broken anyway, so no matter.  Still want to provide Iqxy in case 
     493  // the user model wants full control of orientation/magnetism. 
     494  #if defined(HAVE_PSI) 
     495    const double theta = values[details->theta_par+2]; 
     496    const double phi = values[details->theta_par+3]; 
     497    const double psi = values[details->theta_par+4]; 
     498    double weight; 
     499    #define APPLY_PROJECTION() do { \ 
     500      local_values.table.theta = theta; \ 
     501      local_values.table.phi = phi; \ 
     502      local_values.table.psi = psi; \ 
     503      weight=weight0; \ 
     504    } while (0) 
     505  #elif defined(HAVE_THETA) 
     506    const double theta = values[details->theta_par+2]; 
     507    const double phi = values[details->theta_par+3]; 
     508    double weight; 
     509    #define APPLY_PROJECTION() do { \ 
     510      local_values.table.theta = theta; \ 
     511      local_values.table.phi = phi; \ 
     512      weight=weight0; \ 
     513    } while (0) 
     514  #else 
     515    #define APPLY_PROJECTION() const double weight=weight0 
     516  #endif 
    479517#else // !spherosymmetric projection 
    480518  // Grab the "view" angles (theta, phi, psi) from the initial parameter table. 
  • sasmodels/kernelpy.py

    r2d81cfe r108e70e  
    2626# pylint: enable=unused-import 
    2727 
     28logger = logging.getLogger(__name__) 
     29 
    2830class PyModel(KernelModel): 
    2931    """ 
     
    3133    """ 
    3234    def __init__(self, model_info): 
    33         # Make sure Iq and Iqxy are available and vectorized 
     35        # Make sure Iq is available and vectorized 
    3436        _create_default_functions(model_info) 
    3537        self.info = model_info 
    3638        self.dtype = np.dtype('d') 
    37         logging.info("load python model " + self.info.name) 
     39        logger.info("load python model " + self.info.name) 
    3840 
    3941    def make_kernel(self, q_vectors): 
    4042        q_input = PyInput(q_vectors, dtype=F64) 
    41         kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 
    42         return PyKernel(kernel, self.info, q_input) 
     43        return PyKernel(self.info, q_input) 
    4344 
    4445    def release(self): 
     
    8990    Callable SAS kernel. 
    9091 
    91     *kernel* is the DllKernel object to call. 
     92    *kernel* is the kernel object to call. 
    9293 
    9394    *model_info* is the module information 
     
    104105    Call :meth:`release` when done with the kernel instance. 
    105106    """ 
    106     def __init__(self, kernel, model_info, q_input): 
     107    def __init__(self, model_info, q_input): 
    107108        # type: (callable, ModelInfo, List[np.ndarray]) -> None 
    108109        self.dtype = np.dtype('d') 
     
    110111        self.q_input = q_input 
    111112        self.res = np.empty(q_input.nq, q_input.dtype) 
    112         self.kernel = kernel 
    113113        self.dim = '2d' if q_input.is_2d else '1d' 
    114114 
     
    159159        # Generate a closure which calls the form_volume if it exists. 
    160160        form_volume = model_info.form_volume 
    161         self._volume = ((lambda: form_volume(*volume_args)) if form_volume 
    162                         else (lambda: 1.0)) 
     161        self._volume = ((lambda: form_volume(*volume_args)) if form_volume else 
     162                        (lambda: 1.0)) 
    163163 
    164164    def __call__(self, call_details, values, cutoff, magnetic): 
     
    261261    any functions that are not already marked as vectorized. 
    262262    """ 
     263    # Note: must call create_vector_Iq before create_vector_Iqxy 
    263264    _create_vector_Iq(model_info) 
    264     _create_vector_Iqxy(model_info)  # call create_vector_Iq() first 
     265    _create_vector_Iqxy(model_info) 
    265266 
    266267 
     
    280281        model_info.Iq = vector_Iq 
    281282 
     283 
    282284def _create_vector_Iqxy(model_info): 
    283285    """ 
    284286    Define Iqxy as a vector function if it exists, or default it from Iq(). 
    285287    """ 
    286     Iq, Iqxy = model_info.Iq, model_info.Iqxy 
     288    Iqxy = getattr(model_info, 'Iqxy', None) 
    287289    if callable(Iqxy): 
    288290        if not getattr(Iqxy, 'vectorized', False): 
     
    295297            vector_Iqxy.vectorized = True 
    296298            model_info.Iqxy = vector_Iqxy 
    297     elif callable(Iq): 
     299    else: 
    298300        #print("defaulting Iqxy") 
    299301        # Iq is vectorized because create_vector_Iq was already called. 
     302        Iq = model_info.Iq 
    300303        def default_Iqxy(qx, qy, *args): 
    301304            """ 
  • sasmodels/modelinfo.py

    r2d81cfe r108e70e  
    3737 
    3838# assumptions about common parameters exist throughout the code, such as: 
    39 # (1) kernel functions Iq, Iqxy, form_volume, ... don't see them 
     39# (1) kernel functions Iq, Iqxy, Iqac, Iqabc, form_volume, ... don't see them 
    4040# (2) kernel drivers assume scale is par[0] and background is par[1] 
    4141# (3) mixture models drop the background on components and replace the scale 
     
    256256 
    257257    *type* indicates how the parameter will be used.  "volume" parameters 
    258     will be used in all functions.  "orientation" parameters will be used 
    259     in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in 
    260     *Imagnetic* only.  If *type* is the empty string, the parameter will 
     258    will be used in all functions.  "orientation" parameters are not passed, 
     259    but will be used to convert from *qx*, *qy* to *qa*, *qb*, *qc* in calls to 
     260    *Iqxy* and *Imagnetic*.  If *type* is the empty string, the parameter will 
    261261    be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters 
    262262    can automatically be promoted to magnetic parameters, each of which 
     
    386386      with vector parameter p sent as p[]. 
    387387 
    388     * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 
    389       function, with vector parameter p sent as p[]. 
    390  
    391388    * *form_volume_parameters* is the list of parameters to the form_volume(...) 
    392389      function, with vector parameter p sent as p[]. 
     
    443440        self.iq_parameters = [p for p in self.kernel_parameters 
    444441                              if p.type not in ('orientation', 'magnetic')] 
    445         # note: orientation no longer sent to Iqxy, so its the same as 
    446         #self.iqxy_parameters = [p for p in self.kernel_parameters 
    447         #                        if p.type != 'magnetic'] 
     442        self.orientation_parameters = [p for p in self.kernel_parameters 
     443                                       if p.type == 'orientation'] 
    448444        self.form_volume_parameters = [p for p in self.kernel_parameters 
    449445                                       if p.type == 'volume'] 
     
    490486                if p.type != 'orientation': 
    491487                    raise TypeError("psi must be an orientation parameter") 
     488            elif p.type == 'orientation': 
     489                raise TypeError("only theta, phi and psi can be orientation parameters") 
    492490        if theta >= 0 and phi >= 0: 
     491            last_par = len(self.kernel_parameters) - 1 
    493492            if phi != theta+1: 
    494493                raise TypeError("phi must follow theta") 
    495494            if psi >= 0 and psi != phi+1: 
    496495                raise TypeError("psi must follow phi") 
     496            #if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par): 
     497            #    raise TypeError("orientation parameters must appear at the " 
     498            #                    "end of the parameter table") 
    497499        elif theta >= 0 or phi >= 0 or psi >= 0: 
    498500            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
     
    715717 
    716718 
     719#: Set of variables defined in the model that might contain C code 
     720C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
     721 
    717722def _find_source_lines(model_info, kernel_module): 
    718723    # type: (ModelInfo, ModuleType) -> None 
     
    720725    Identify the location of the C source inside the model definition file. 
    721726 
    722     This code runs through the source of the kernel module looking for 
    723     lines that start with 'Iq', 'Iqxy' or 'form_volume'.  Clearly there are 
    724     all sorts of reasons why this might not work (e.g., code commented out 
    725     in a triple-quoted line block, code built using string concatenation, 
    726     or code defined in the branch of an 'if' block), but it should work 
    727     properly in the 95% case, and getting the incorrect line number will 
    728     be harmless. 
    729     """ 
    730     # Check if we need line numbers at all 
    731     if callable(model_info.Iq): 
    732         return None 
    733  
    734     if (model_info.Iq is None 
    735             and model_info.Iqxy is None 
    736             and model_info.Imagnetic is None 
    737             and model_info.form_volume is None): 
     727    This code runs through the source of the kernel module looking for lines 
     728    that contain C code (because it is a c function definition).  Clearly 
     729    there are all sorts of reasons why this might not work (e.g., code 
     730    commented out in a triple-quoted line block, code built using string 
     731    concatenation, code defined in the branch of an 'if' block, code imported 
     732    from another file), but it should work properly in the 95% case, and for 
     733    the remainder, getting the incorrect line number will merely be 
     734    inconvenient. 
     735    """ 
     736    # Only need line numbers if we are creating a C module and the C symbols 
     737    # are defined. 
     738    if (callable(model_info.Iq) 
     739            or not any(hasattr(model_info, s) for s in C_SYMBOLS)): 
    738740        return 
    739741 
    740     # find the defintion lines for the different code blocks 
     742    # load the module source if we can 
    741743    try: 
    742744        source = inspect.getsource(kernel_module) 
    743745    except IOError: 
    744746        return 
    745     for k, v in enumerate(source.split('\n')): 
    746         if v.startswith('Imagnetic'): 
    747             model_info._Imagnetic_line = k+1 
    748         elif v.startswith('Iqxy'): 
    749             model_info._Iqxy_line = k+1 
    750         elif v.startswith('Iq'): 
    751             model_info._Iq_line = k+1 
    752         elif v.startswith('form_volume'): 
    753             model_info._form_volume_line = k+1 
    754  
     747 
     748    # look for symbol at the start of the line 
     749    for lineno, line in enumerate(source.split('\n')): 
     750        for name in C_SYMBOLS: 
     751            if line.startswith(name): 
     752                # Add 1 since some compilers complain about "#line 0" 
     753                model_info.lineno[name] = lineno + 1 
     754                break 
    755755 
    756756def make_model_info(kernel_module): 
     
    761761    Fill in default values for parts of the module that are not provided. 
    762762 
    763     Note: vectorized Iq and Iqxy functions will be created for python 
     763    Note: vectorized Iq and Iqac/Iqabc functions will be created for python 
    764764    models when the model is first called, not when the model is loaded. 
    765765    """ 
     
    790790    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    791791    info.source = getattr(kernel_module, 'source', []) 
     792    info.c_code = getattr(kernel_module, 'c_code', None) 
    792793    # TODO: check the structure of the tests 
    793794    info.tests = getattr(kernel_module, 'tests', []) 
     
    797798    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    798799    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 
     800    info.Iqac = getattr(kernel_module, 'Iqac', None) # type: ignore 
     801    info.Iqabc = getattr(kernel_module, 'Iqabc', None) # type: ignore 
    799802    info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore 
    800803    info.profile = getattr(kernel_module, 'profile', None) # type: ignore 
     
    811814    info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 
    812815 
     816    if callable(info.Iq) and parameters.has_2d: 
     817        raise ValueError("oriented python models not supported") 
     818 
     819    info.lineno = {} 
    813820    _find_source_lines(info, kernel_module) 
    814821 
     
    824831 
    825832    The structure should be mostly static, other than the delayed definition 
    826     of *Iq* and *Iqxy* if they need to be defined. 
     833    of *Iq*, *Iqac* and *Iqabc* if they need to be defined. 
    827834    """ 
    828835    #: Full path to the file defining the kernel, if any. 
     
    906913    structure_factor = None # type: bool 
    907914    #: List of C source files used to define the model.  The source files 
    908     #: should define the *Iq* function, and possibly *Iqxy*, though a default 
    909     #: *Iqxy = Iq(sqrt(qx**2+qy**2)* will be created if no *Iqxy* is provided. 
    910     #: Files containing the most basic functions must appear first in the list, 
    911     #: followed by the files that use those functions.  Form factors are 
    912     #: indicated by providing a :attr:`ER` function. 
     915    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 
     916    #: model defines orientation parameters. Files containing the most basic 
     917    #: functions must appear first in the list, followed by the files that 
     918    #: use those functions.  Form factors are indicated by providing 
     919    #: an :attr:`ER` function. 
    913920    source = None           # type: List[str] 
    914921    #: The set of tests that must pass.  The format of the tests is described 
     
    935942    #: See :attr:`ER` for details on the parameters. 
    936943    VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]] 
     944    #: Arbitrary C code containing supporting functions, etc., to be inserted 
     945    #: after everything in source.  This can include Iq and Iqxy functions with 
     946    #: the full function signature, including all parameters. 
     947    c_code = None 
    937948    #: Returns the form volume for python-based models.  Form volume is needed 
    938949    #: for volume normalization in the polydispersity integral.  If no 
     
    955966    #: include the decimal point. See :mod:`generate` for more details. 
    956967    Iq = None               # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    957     #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
    958     Iqxy = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     968    #: Returns *I(qab, qc, a, b, ...)*.  The interface follows :attr:`Iq`. 
     969    Iqac = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     970    #: Returns *I(qa, qb, qc, a, b, ...)*.  The interface follows :attr:`Iq`. 
     971    Iqabc = None            # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    959972    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
    960973    Imagnetic = None        # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     
    972985    #: Returns a random parameter set for the model 
    973986    random = None           # type: Optional[Callable[[], Dict[str, float]]] 
    974  
    975     # line numbers within the python file for bits of C source, if defined 
    976     # NB: some compilers fail with a "#line 0" directive, so default to 1. 
    977     _Imagnetic_line = 1 
    978     _Iqxy_line = 1 
    979     _Iq_line = 1 
    980     _form_volume_line = 1 
    981  
     987    #: Line numbers for symbols defining C code 
     988    lineno = None           # type: Dict[str, int] 
    982989 
    983990    def __init__(self): 
  • sasmodels/models/_spherepy.py

    ref07e95 r108e70e  
    8888Iq.vectorized = True  # Iq accepts an array of q values 
    8989 
    90 def Iqxy(qx, qy, sld, sld_solvent, radius): 
    91     return Iq(sqrt(qx ** 2 + qy ** 2), sld, sld_solvent, radius) 
    92 Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values 
    93  
    9490def sesans(z, sld, sld_solvent, radius): 
    9591    """ 
  • sasmodels/models/barbell.c

    r74768cb r108e70e  
    9090 
    9191static double 
    92 Iqxy(double qab, double qc, 
     92Iqac(double qab, double qc, 
    9393    double sld, double solvent_sld, 
    9494    double radius_bell, double radius, double length) 
  • sasmodels/models/bcc_paracrystal.c

    r74768cb r108e70e  
    107107 
    108108 
    109 static double Iqxy(double qa, double qb, double qc, 
     109static double Iqabc(double qa, double qb, double qc, 
    110110    double dnn, double d_factor, double radius, 
    111111    double sld, double solvent_sld) 
  • sasmodels/models/capped_cylinder.c

    r74768cb r108e70e  
    115115 
    116116static double 
    117 Iqxy(double qab, double qc, 
     117Iqac(double qab, double qc, 
    118118    double sld, double solvent_sld, double radius, 
    119119    double radius_cap, double length) 
  • sasmodels/models/core_shell_bicelle.c

    r74768cb r108e70e  
    6767 
    6868static double 
    69 Iqxy(double qab, double qc, 
     69Iqac(double qab, double qc, 
    7070    double radius, 
    7171    double thick_rim, 
  • sasmodels/models/core_shell_bicelle_elliptical.c

    rd4db147 r108e70e  
    7171 
    7272static double 
    73 Iqxy(double qa, double qb, double qc, 
     73Iqabc(double qa, double qb, double qc, 
    7474    double r_minor, 
    7575    double x_core, 
  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    rd4db147 r108e70e  
    7979 
    8080static double 
    81 Iqxy(double qa, double qb, double qc, 
     81Iqabc(double qa, double qb, double qc, 
    8282          double r_minor, 
    8383          double x_core, 
     
    114114    return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma)); 
    115115} 
    116  
  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    r110f69c r108e70e  
    149149    ["sld_rim",        "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Cylinder rim scattering length density"], 
    150150    ["sld_solvent",    "1e-6/Ang^2", 6, [-inf, inf], "sld",         "Solvent scattering length density"], 
     151    ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"], 
    151152    ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 
    152153    ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"], 
    153154    ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about cylinder axis"], 
    154     ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"] 
    155155    ] 
    156156 
  • sasmodels/models/core_shell_cylinder.c

    r74768cb r108e70e  
    4848 
    4949 
    50 double Iqxy(double qab, double qc, 
     50double Iqac(double qab, double qc, 
    5151    double core_sld, 
    5252    double shell_sld, 
  • sasmodels/models/core_shell_ellipsoid.c

    r74768cb r108e70e  
    7575 
    7676static double 
    77 Iqxy(double qab, double qc, 
     77Iqac(double qab, double qc, 
    7878    double radius_equat_core, 
    7979    double x_core, 
  • sasmodels/models/core_shell_parallelepiped.c

    ra261a83 r108e70e  
    102102 
    103103static double 
    104 Iqxy(double qa, double qb, double qc, 
     104Iqabc(double qa, double qb, double qc, 
    105105    double core_sld, 
    106106    double arim_sld, 
  • sasmodels/models/cylinder.c

    r74768cb r108e70e  
    4545 
    4646static double 
    47 Iqxy(double qab, double qc, 
     47Iqac(double qab, double qc, 
    4848    double sld, 
    4949    double solvent_sld, 
  • sasmodels/models/ellipsoid.c

    r74768cb r108e70e  
    3939 
    4040static double 
    41 Iqxy(double qab, double qc, 
     41Iqac(double qab, double qc, 
    4242    double sld, 
    4343    double sld_solvent, 
  • sasmodels/models/elliptical_cylinder.c

    rd4db147 r108e70e  
    5454 
    5555static double 
    56 Iqxy(double qa, double qb, double qc, 
     56Iqabc(double qa, double qb, double qc, 
    5757     double radius_minor, double r_ratio, double length, 
    5858     double sld, double solvent_sld) 
  • sasmodels/models/fcc_paracrystal.c

    r74768cb r108e70e  
    8080 
    8181 
    82 static double Iqxy(double qa, double qb, double qc, 
     82static double Iqabc(double qa, double qb, double qc, 
    8383    double dnn, double d_factor, double radius, 
    8484    double sld, double solvent_sld) 
  • sasmodels/models/hollow_cylinder.c

    r74768cb r108e70e  
    5252 
    5353static double 
    54 Iqxy(double qab, double qc, 
     54Iqac(double qab, double qc, 
    5555    double radius, double thickness, double length, 
    5656    double sld, double solvent_sld) 
  • sasmodels/models/hollow_rectangular_prism.c

    r74768cb r108e70e  
    8585} 
    8686 
    87 double Iqxy(double qa, double qb, double qc, 
     87double Iqabc(double qa, double qb, double qc, 
    8888    double sld, 
    8989    double solvent_sld, 
  • sasmodels/models/line.py

    r2d81cfe r108e70e  
    5757Iq.vectorized = True # Iq accepts an array of q values 
    5858 
     59 
    5960def Iqxy(qx, qy, *args): 
    6061    """ 
     
    6970 
    7071Iqxy.vectorized = True  # Iqxy accepts an array of qx qy values 
     72 
     73# uncomment the following to test Iqxy in C models 
     74#del Iq, Iqxy 
     75#c_code = """ 
     76#static double Iq(double q, double b, double m) { return m*q+b; } 
     77#static double Iqxy(double qx, double qy, double b, double m) 
     78#{ return (m*qx+b)*(m*qy+b); } 
     79#""" 
    7180 
    7281def random(): 
  • sasmodels/models/parallelepiped.c

    r74768cb r108e70e  
    5353 
    5454static double 
    55 Iqxy(double qa, double qb, double qc, 
     55Iqabc(double qa, double qb, double qc, 
    5656    double sld, 
    5757    double solvent_sld, 
  • sasmodels/models/rectangular_prism.c

    r74768cb r108e70e  
    7171 
    7272 
    73 double Iqxy(double qa, double qb, double qc, 
     73double Iqabc(double qa, double qb, double qc, 
    7474    double sld, 
    7575    double solvent_sld, 
  • sasmodels/models/sc_paracrystal.c

    r74768cb r108e70e  
    8282 
    8383static double 
    84 Iqxy(double qa, double qb, double qc, 
     84Iqabc(double qa, double qb, double qc, 
    8585    double dnn, double d_factor, double radius, 
    8686    double sld, double solvent_sld) 
  • sasmodels/models/stacked_disks.c

    r74768cb r108e70e  
    142142 
    143143static double 
    144 Iqxy(double qab, double qc, 
     144Iqac(double qab, double qc, 
    145145    double thick_core, 
    146146    double thick_layer, 
  • sasmodels/models/triaxial_ellipsoid.c

    r74768cb r108e70e  
    4646 
    4747static double 
    48 Iqxy(double qa, double qb, double qc, 
     48Iqabc(double qa, double qb, double qc, 
    4949    double sld, 
    5050    double sld_solvent, 
Note: See TracChangeset for help on using the changeset viewer.