Changeset 8698a0d in sasmodels


Ignore:
Timestamp:
Oct 17, 2017 11:53:01 PM (7 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:
32f87a5
Parents:
becded3
git-author:
Paul Kienzle <pkienzle@…> (10/17/17 18:23:09)
git-committer:
Paul Kienzle <pkienzle@…> (10/17/17 23:53:01)
Message:

revise api for oriented shapes, allowing jitter in the frame of the sample

Location:
sasmodels
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    rb76191e r8698a0d  
    9191 
    9292    === precision options === 
    93     -calc=default uses the default calcution precision 
     93    -engine=default uses the default calcution precision 
    9494    -single/-double/-half/-fast sets an OpenCL calculation engine 
    9595    -single!/-double!/-quad! sets an OpenMP calculation engine 
    … …  
    110110vary from 64-bit to 128-bit, with 80-bit floats being common (1e-19 precision). 
    111111On unix and mac you may need single quotes around the DLL computation 
    112 engines, such as -calc='single!,double!' since !, is treated as a history 
     112engines, such as -engine='single!,double!' since !, is treated as a history 
    113113expansion request in the shell. 
    114114 
    … …  
    122122 
    123123    # compare single and double precision calculation for a barbell 
    124     sascomp barbell -calc=single,double 
     124    sascomp barbell -engine=single,double 
    125125 
    126126    # generate 10 random lorentz models, with seed=27 
    … …  
    131131 
    132132    # model timing test requires multiple evals to perform the estimate 
    133     sascomp pringle -calc=single,double -timing=100,100 -noplot 
     133    sascomp pringle -engine=single,double -timing=100,100 -noplot 
    134134""" 
    135135 
    … …  
    10171017 
    10181018    # Precision options 
    1019     'calc=', 
     1019    'engine=', 
    10201020    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!', 
    10211021    'sasview',  # TODO: remove sasview 3.x support 
    … …  
    11671167        elif arg.startswith('-title='):    opts['title'] = arg[7:] 
    11681168        elif arg.startswith('-data='):     opts['datafile'] = arg[6:] 
    1169         elif arg.startswith('-calc='):     opts['engine'] = arg[6:] 
     1169        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
    11701170        elif arg.startswith('-neval='):    opts['evals'] = arg[7:] 
    11711171        elif arg == '-random':  opts['seed'] = np.random.randint(1000000) 

  • sasmodels/details.py

    rf39759c r8698a0d  
    219219 
    220220ZEROS = tuple([0.]*31) 
    221 def make_kernel_args(kernel, pairs): 
     221def make_kernel_args(kernel, mesh): 
    222222    # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool] 
    223223    """ 
    224     Converts (value, weight) pairs into parameters for the kernel call. 
     224    Converts (value, dispersity, weight) for each parameter into kernel pars. 
    225225 
    226226    Returns a CallDetails object indicating the polydispersity, a data object 
    … …  
    231231    npars = kernel.info.parameters.npars 
    232232    nvalues = kernel.info.parameters.nvalues 
    233     scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs] 
     233    scalars = [value for value, dispersity, weight in mesh] 
    234234    # skipping scale and background when building values and weights 
    235     values, weights = zip(*pairs[2:npars+2]) if npars else ((), ()) 
    236     weights = correct_theta_weights(kernel.info.parameters, values, weights) 
     235    values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 
     236    #weights = correct_theta_weights(kernel.info.parameters, values, weights) 
    237237    length = np.array([len(w) for w in weights]) 
    238238    offset = np.cumsum(np.hstack((0, length))) 
    239239    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
    240240    # Pad value array to a 32 value boundaryd 
    241     data_len = nvalues + 2*sum(len(v) for v in values) 
     241    data_len = nvalues + 2*sum(len(v) for v in dispersity) 
    242242    extra = (32 - data_len%32)%32 
    243     data = np.hstack((scalars,) + values + weights + ZEROS[:extra]) 
     243    data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra]) 
    244244    data = data.astype(kernel.dtype) 
    245245    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
    … …  
    294294 
    295295 
    296 def dispersion_mesh(model_info, pars): 
     296def dispersion_mesh(model_info, mesh): 
    297297    # type: (ModelInfo) -> Tuple[List[np.ndarray], List[np.ndarray]] 
    298298    """ 
    299299    Create a mesh grid of dispersion parameters and weights. 
    300300 
    301     pars is a list of pairs (values, weights), where the values are the 
    302     individual parameter values at which to evaluate the polydispersity 
    303     distribution and weights are the weights associated with each value. 
     301    *mesh* is a list of (value, dispersity, weights), where the values 
     302    are the individual parameter values, and (dispersity, weights) is 
     303    the distribution of parameter values. 
    304304 
    305305    Only the volume parameters should be included in this list.  Orientation 
    306306    parameters do not affect the calculation of effective radius or volume 
    307     ratio. 
     307    ratio.  This is convenient since it avoids the distinction between 
     308    value and dispersity that is present in orientation parameters but not 
     309    shape parameters. 
    308310 
    309311    Returns [p1,p2,...],w where pj is a vector of values for parameter j 
    … …  
    311313    parameter set in the vector. 
    312314    """ 
    313     value, weight = zip(*pars) 
     315    value, dispersity, weight = zip(*mesh) 
    314316    #weight = [w if len(w)>0 else [1.] for w in weight] 
    315317    weight = np.vstack([v.flatten() for v in meshgrid(*weight)]) 
    316318    weight = np.prod(weight, axis=0) 
    317     value = [v.flatten() for v in meshgrid(*value)] 
     319    dispersity = [v.flatten() for v in meshgrid(*dispersity)] 
    318320    lengths = [par.length for par in model_info.parameters.kernel_parameters 
    319321               if par.type == 'volume'] 
    … …  
    322324        offset = 0 
    323325        for n in lengths: 
    324             pars.append(np.vstack(value[offset:offset+n]) 
    325                         if n > 1 else value[offset]) 
     326            pars.append(np.vstack(dispersity[offset:offset+n]) 
     327                        if n > 1 else dispersity[offset]) 
    326328            offset += n 
    327         value = pars 
    328     return value, weight 
     329        dispersity = pars 
     330    return dispersity, weight 

  • sasmodels/direct_model.py

    rd1ff3a5 r8698a0d  
    6666 
    6767    #print("pars",[p.id for p in parameters.call_parameters]) 
    68     vw_pairs = [(get_weights(p, pars) if active(p.name) 
    69                  else ([pars.get(p.name, p.default)], [1.0])) 
    70                 for p in parameters.call_parameters] 
    71  
    72     call_details, values, is_magnetic = make_kernel_args(calculator, vw_pairs) 
     68    mesh = [get_weights(p, pars, active(p.name)) 
     69            for p in parameters.call_parameters] 
     70 
     71    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
    7372    #print("values:", values) 
    7473    return calculator(call_details, values, cutoff, is_magnetic) 
    … …  
    130129 
    131130 
    132 def get_weights(parameter, values): 
     131def get_weights(parameter, values, active=True): 
    133132    # type: (Parameter, Dict[str, float]) -> Tuple[np.ndarray, np.ndarray] 
    134133    """ 
    … …  
    140139    """ 
    141140    value = float(values.get(parameter.name, parameter.default)) 
    142     relative = parameter.relative_pd 
    143     limits = parameter.limits 
    144     disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
    145141    npts = values.get(parameter.name+'_pd_n', 0) 
    146142    width = values.get(parameter.name+'_pd', 0.0) 
    147     nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
    148     if npts == 0 or width == 0: 
    149         return [value], [1.0] 
    150     value, weight = weights.get_weights( 
    151         disperser, npts, width, nsigma, value, limits, relative) 
    152     return value, weight / np.sum(weight) 
     143    if npts == 0 or width == 0 or not active: 
     144        # Note: orientation parameters have the viewing angle as the parameter 
     145        # value and the jitter in the distribution, so be sure to set the 
     146        # empty pd for orientation parameters to 0. 
     147        pd = [value if parameter.type != 'orientation' else 0.0], [1.0] 
     148    else: 
     149        relative = parameter.relative_pd 
     150        limits = parameter.limits 
     151        disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
     152        nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
     153        pd = weights.get_weights(disperser, npts, width, nsigma, 
     154                                value, limits, relative) 
     155    return value, pd[0], pd[1] 
    153156 
    154157 

  • sasmodels/generate.py

    r30b60d2 r8698a0d  
    714714    # Define the parameter table 
    715715    # TODO: plug in current line number 
    716     source.append('#line 542 "sasmodels/generate.py"') 
     716    source.append('#line 716 "sasmodels/generate.py"') 
    717717    source.append("#define PARAMETER_TABLE \\") 
    718718    source.append("\\\n".join(p.as_definition() 
    … …  
    730730    source.append(call_volume) 
    731731 
    732     refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 
    733     call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 
     732    model_refs = _call_pars("_v.", partable.iq_parameters) 
     733    pars = ",".join(["_q"] + model_refs) 
     734    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
    734735    if _have_Iqxy(user_code) or isinstance(model_info.Iqxy, str): 
    735         # Call 2D model 
    736         refs = ["_q[2*_i]", "_q[2*_i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 
    737         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 
     736        if partable.is_asymmetric: 
     737            pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
     738            call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
     739            clear_iqxy = "#undef CALL_IQ_ABC" 
     740        else: 
     741            pars = ",".join(["_qa", "_qc"] + model_refs) 
     742            call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 
     743            clear_iqxy = "#undef CALL_IQ_AC" 
    738744    else: 
    739         # Call 1D model with sqrt(qx^2 + qy^2) 
    740         #warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
    741         # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
    742         pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:] 
    743         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
     745        pars = ",".join(["_qa"] + model_refs) 
     746        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
     747        clear_iqxy = "#undef CALL_IQ_A" 
    744748 
    745749    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    … …  
    757761    # TODO: allow mixed python/opencl kernels? 
    758762 
    759     ocl = kernels(ocl_code, call_iq, call_iqxy, model_info.name) 
    760     dll = kernels(dll_code, call_iq, call_iqxy, model_info.name) 
     763    ocl = kernels(ocl_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
     764    dll = kernels(dll_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    761765    result = { 
    762766        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    … …  
    767771 
    768772 
    769 def kernels(kernel, call_iq, call_iqxy, name): 
     773def kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
    770774    # type: ([str,str], str, str, str) -> List[str] 
    771775    code = kernel[0] 
    … …  
    787791        '#line 1 "%s Iqxy"' % path, 
    788792        code, 
    789         "#undef CALL_IQ", 
     793        clear_iqxy, 
    790794        "#undef KERNEL_NAME", 
    791795    ] 
    … …  
    798802        '#line 1 "%s Imagnetic"' % path, 
    799803        code, 
     804        clear_iqxy, 
    800805        "#undef MAGNETIC", 
    801         "#undef CALL_IQ", 
    802806        "#undef KERNEL_NAME", 
    803807    ] 

  • sasmodels/kernel_header.c

    r73cbc5b r8698a0d  
    8787 
    8888#if defined(NEED_EXPM1) 
     89   // TODO: precision is a half digit lower than numpy on mac in [1e-7, 0.5] 
     90   // Run "explore/precision.py sas_expm1" to see this (may have to fiddle 
     91   // the xrange for log to see the complete range). 
    8992   static SAS_DOUBLE expm1(SAS_DOUBLE x_in) { 
    9093      double x = (double)x_in;  // go back to float for single precision kernels 
    … …  
    147150inline double cube(double x) { return x*x*x; } 
    148151inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    149  
    150 // To rotate from the canonical position to theta, phi, psi, first rotate by 
    151 // psi about the major axis, oriented along z, which is a rotation in the 
    152 // detector plane xy. Next rotate by theta about the y axis, aligning the major 
    153 // axis in the xz plane. Finally, rotate by phi in the detector plane xy. 
    154 // To compute the scattering, undo these rotations in reverse order: 
    155 //     rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi 
    156 // The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit 
    157 // vector in the q direction. 
    158 // To change between counterclockwise and clockwise rotation, change the 
    159 // sign of phi and psi. 
    160  
    161 #if 1 
    162 //think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
    163 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
    164     SINCOS(phi*M_PI_180, sn, cn); \ 
    165     q = sqrt(qx*qx + qy*qy); \ 
    166     cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
    167     sn = sqrt(1 - cn*cn); \ 
    168     } while (0) 
    169 #else 
    170 // SasView 3.x definition of orientation 
    171 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
    172     SINCOS(theta*M_PI_180, sn, cn); \ 
    173     q = sqrt(qx*qx + qy*qy);\ 
    174     cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \ 
    175     sn = sqrt(1 - cn*cn); \ 
    176     } while (0) 
    177 #endif 
    178  
    179 #if 1 
    180 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \ 
    181     q = sqrt(qx*qx + qy*qy); \ 
    182     const double qxhat = qx/q; \ 
    183     const double qyhat = qy/q; \ 
    184     double sin_theta, cos_theta; \ 
    185     double sin_phi, cos_phi; \ 
    186     double sin_psi, cos_psi; \ 
    187     SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
    188     SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
    189     SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
    190     xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \ 
    191          + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \ 
    192     yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \ 
    193          + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \ 
    194     zhat = qxhat*(-sin_theta*cos_phi) \ 
    195          + qyhat*(-sin_theta*sin_phi); \ 
    196     } while (0) 
    197 #else 
    198 // SasView 3.x definition of orientation 
    199 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
    200     q = sqrt(qx*qx + qy*qy); \ 
    201     const double qxhat = qx/q; \ 
    202     const double qyhat = qy/q; \ 
    203     double sin_theta, cos_theta; \ 
    204     double sin_phi, cos_phi; \ 
    205     double sin_psi, cos_psi; \ 
    206     SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
    207     SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
    208     SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
    209     cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ 
    210     cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ 
    211     cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
    212     } while (0) 
    213 #endif 

  • sasmodels/kernel_iq.c

    rd4c33d6 r8698a0d  
    2525    int32_t num_weights;        // total length of the weights vector 
    2626    int32_t num_active;         // number of non-trivial pd loops 
    27     int32_t theta_par;          // id of spherical correction variable (not used) 
     27    int32_t theta_par;          // id of first orientation variable 
    2828} ProblemDetails; 
    2929 
    … …  
    6969    return sld + perp*p; 
    7070} 
    71  
    72 #endif // MAGNETIC 
     71//#endif // MAGNETIC 
     72 
     73// TODO: way too hackish 
     74// For the 1D kernel, CALL_IQ_[A,AC,ABC] and MAGNETIC are not defined 
     75// so view_direct *IS NOT* included 
     76// For the 2D kernel, CALL_IQ_[A,AC,ABC] is defined but MAGNETIC is not 
     77// so view_direct *IS* included 
     78// For the magnetic kernel, CALL_IQ_[A,AC,ABC] is defined, but so is MAGNETIC 
     79// so view_direct *IS NOT* included 
     80#else // !MAGNETIC 
     81 
     82// ===== Implement jitter in orientation ===== 
     83// To change the definition of the angles, run explore/angles.py, which 
     84// uses sympy to generate the equations. 
     85 
     86#if defined(CALL_IQ_AC) // oriented symmetric 
     87static double 
     88view_direct(double qx, double qy, 
     89             double theta, double phi, 
     90             ParameterTable table) 
     91{ 
     92    double sin_theta, cos_theta; 
     93    double sin_phi, cos_phi; 
     94 
     95    // reverse view 
     96    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     97    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     98    const double qa = qx*cos_phi*cos_theta + qy*sin_phi*cos_theta; 
     99    const double qb = -qx*sin_phi + qy*cos_phi; 
     100    const double qc = qx*sin_theta*cos_phi + qy*sin_phi*sin_theta; 
     101 
     102    // reverse jitter after view 
     103    SINCOS(table.theta*M_PI_180, sin_theta, cos_theta); 
     104    SINCOS(table.phi*M_PI_180, sin_phi, cos_phi); 
     105    const double dqc = qa*sin_theta - qb*sin_phi*cos_theta + qc*cos_phi*cos_theta; 
     106 
     107    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
     108    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
     109 
     110    return CALL_IQ_AC(dqa, dqc, table); 
     111} 
     112 
     113#elif defined(CALL_IQ_ABC) // oriented asymmetric 
     114 
     115static double 
     116view_direct(double qx, double qy, 
     117             double theta, double phi, double psi, 
     118             ParameterTable table) 
     119{ 
     120    double sin_theta, cos_theta; 
     121    double sin_phi, cos_phi; 
     122    double sin_psi, cos_psi; 
     123 
     124    // reverse view 
     125    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     126    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     127    SINCOS(psi*M_PI_180, sin_psi, cos_psi); 
     128    const double qa = qx*(sin_phi*sin_psi + cos_phi*cos_psi*cos_theta) + qy*(sin_phi*cos_psi*cos_theta - sin_psi*cos_phi); 
     129    const double qb = qx*(-sin_phi*cos_psi + sin_psi*cos_phi*cos_theta) + qy*(sin_phi*sin_psi*cos_theta + cos_phi*cos_psi); 
     130    const double qc = qx*sin_theta*cos_phi + qy*sin_phi*sin_theta; 
     131 
     132    // reverse jitter after view 
     133    SINCOS(table.theta*M_PI_180, sin_theta, cos_theta); 
     134    SINCOS(table.phi*M_PI_180, sin_phi, cos_phi); 
     135    SINCOS(table.psi*M_PI_180, sin_psi, cos_psi); 
     136    const double dqa = qa*cos_psi*cos_theta + qb*(sin_phi*sin_theta*cos_psi - sin_psi*cos_phi) + qc*(-sin_phi*sin_psi - sin_theta*cos_phi*cos_psi); 
     137    const double dqb = qa*sin_psi*cos_theta + qb*(sin_phi*sin_psi*sin_theta + cos_phi*cos_psi) + qc*(sin_phi*cos_psi - sin_psi*sin_theta*cos_phi); 
     138    const double dqc = qa*sin_theta - qb*sin_phi*cos_theta + qc*cos_phi*cos_theta; 
     139 
     140    return CALL_IQ_ABC(dqa, dqb, dqc, table); 
     141} 
     142/* TODO: use precalculated jitter for faster 2D calcs on DLL. 
     143static void 
     144view_precalc( 
     145    double theta, double phi, double psi, 
     146    ParameterTable table, 
     147    double *R11, double *R12, double *R21, 
     148    double *R22, double *R31, double *R32) 
     149{ 
     150    double sin_theta, cos_theta; 
     151    double sin_phi, cos_phi; 
     152    double sin_psi, cos_psi; 
     153 
     154    // reverse view matrix 
     155    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     156    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     157    SINCOS(psi*M_PI_180, sin_psi, cos_psi); 
     158    const double V11 = sin_phi*sin_psi + cos_phi*cos_psi*cos_theta; 
     159    const double V12 = sin_phi*cos_psi*cos_theta - sin_psi*cos_phi; 
     160    const double V21 = -sin_phi*cos_psi + sin_psi*cos_phi*cos_theta; 
     161    const double V22 = sin_phi*sin_psi*cos_theta + cos_phi*cos_psi; 
     162    const double V31 = sin_theta*cos_phi; 
     163    const double V32 = sin_phi*sin_theta; 
     164 
     165    // reverse jitter matrix 
     166    SINCOS(table.theta*M_PI_180, sin_theta, cos_theta); 
     167    SINCOS(table.phi*M_PI_180, sin_phi, cos_phi); 
     168    SINCOS(table.psi*M_PI_180, sin_psi, cos_psi); 
     169    const double J11 = cos_psi*cos_theta; 
     170    const double J12 = sin_phi*sin_theta*cos_psi - sin_psi*cos_phi; 
     171    const double J13 = -sin_phi*sin_psi - sin_theta*cos_phi*cos_psi; 
     172    const double J21 = sin_psi*cos_theta; 
     173    const double J22 = sin_phi*sin_psi*sin_theta + cos_phi*cos_psi; 
     174    const double J23 = sin_phi*cos_psi - sin_psi*sin_theta*cos_phi; 
     175    const double J31 = sin_theta; 
     176    const double J32 = -sin_phi*cos_theta; 
     177    const double J33 = cos_phi*cos_theta; 
     178 
     179    // reverse matrix 
     180    *R11 = J11*V11 + J12*V21 + J13*V31; 
     181    *R12 = J11*V12 + J12*V22 + J13*V32; 
     182    *R21 = J21*V11 + J22*V21 + J23*V31; 
     183    *R22 = J21*V12 + J22*V22 + J23*V32; 
     184    *R31 = J31*V11 + J32*V21 + J33*V31; 
     185    *R32 = J31*V12 + J32*V22 + J33*V32; 
     186} 
     187 
     188static double 
     189view_apply(double qx, double qy, 
     190    double R11, double R12, double R21, double R22, double R31, double R32, 
     191    ParameterTable table) 
     192{ 
     193    const double dqa = R11*qx + R12*qy; 
     194    const double dqb = R21*qx + R22*qy; 
     195    const double dqc = R31*qx + R32*qy; 
     196 
     197    CALL_IQ_ABC(dqa, dqb, dqc, table); 
     198} 
     199*/ 
     200#endif 
     201 
     202#endif // !MAGNETIC 
    73203 
    74204kernel 
    … …  
    105235  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
    106236#endif // MAGNETIC 
     237 
     238#if defined(CALL_IQ_AC) // oriented symmetric 
     239  const double theta = values[details->theta_par+2]; 
     240  const double phi = values[details->theta_par+3]; 
     241#elif defined(CALL_IQ_ABC) // oriented asymmetric 
     242  const double theta = values[details->theta_par+2]; 
     243  const double phi = values[details->theta_par+3]; 
     244  const double psi = values[details->theta_par+4]; 
     245#endif 
    107246 
    108247  // Fill in the initial variables 
    … …  
    275414                  } 
    276415                  #endif 
    277                   scattering += CALL_IQ(q, q_index, local_values.table); 
     416#  if defined(CALL_IQ_A) // unoriented 
     417                  scattering += CALL_IQ_A(sqrt(qsq), local_values.table); 
     418#  elif defined(CALL_IQ_AC) // oriented symmetric 
     419                  scattering += view_direct(qx, qy, theta, phi, local_values.table); 
     420#  elif defined(CALL_IQ_ABC) // oriented asymmetric 
     421                  scattering += view_direct(qx, qy, theta, phi, psi, local_values.table); 
     422#  endif 
    278423                } 
    279424              } 
    280425            } 
    281426          } 
    282 #else  // !MAGNETIC 
    283           const double scattering = CALL_IQ(q, q_index, local_values.table); 
     427#elif defined(CALL_IQ) // 1d, not MAGNETIC 
     428          const double scattering = CALL_IQ(q[q_index], local_values.table); 
     429#else  // 2d data, not magnetic 
     430          const double qx = q[2*q_index]; 
     431          const double qy = q[2*q_index+1]; 
     432#  if defined(CALL_IQ_A) // unoriented 
     433          const double scattering = CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table); 
     434#  elif defined(CALL_IQ_AC) // oriented symmetric 
     435          const double scattering = view_direct(qx, qy, theta, phi, local_values.table); 
     436#  elif defined(CALL_IQ_ABC) // oriented asymmetric 
     437          const double scattering = view_direct(qx, qy, theta, phi, psi, local_values.table); 
     438#  endif 
    284439#endif // !MAGNETIC 
    285440//printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0); 

  • sasmodels/kernel_iq.cl

    rd4c33d6 r8698a0d  
    7070} 
    7171 
    72 #endif // MAGNETIC 
     72//#endif // MAGNETIC 
     73 
     74// TODO: way too hackish 
     75// For the 1D kernel, CALL_IQ_[A,AC,ABC] and MAGNETIC are not defined 
     76// so view_direct *IS NOT* included 
     77// For the 2D kernel, CALL_IQ_[A,AC,ABC] is defined but MAGNETIC is not 
     78// so view_direct *IS* included 
     79// For the magnetic kernel, CALL_IQ_[A,AC,ABC] is defined, but so is MAGNETIC 
     80// so view_direct *IS NOT* included 
     81#else // !MAGNETIC 
     82 
     83// ===== Implement jitter in orientation ===== 
     84// To change the definition of the angles, run explore/angles.py, which 
     85// uses sympy to generate the equations. 
     86 
     87#if defined(CALL_IQ_AC) // oriented symmetric 
     88static double 
     89view_direct(double qx, double qy, 
     90             double theta, double phi, 
     91             ParameterTable table) 
     92{ 
     93    double sin_theta, cos_theta; 
     94    double sin_phi, cos_phi; 
     95 
     96    // reverse view 
     97    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     98    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     99    const double qa = qx*cos_phi*cos_theta + qy*sin_phi*cos_theta; 
     100    const double qb = -qx*sin_phi + qy*cos_phi; 
     101    const double qc = qx*sin_theta*cos_phi + qy*sin_phi*sin_theta; 
     102 
     103    // reverse jitter after view 
     104    SINCOS(table.theta*M_PI_180, sin_theta, cos_theta); 
     105    SINCOS(table.phi*M_PI_180, sin_phi, cos_phi); 
     106    const double dqc = qa*sin_theta - qb*sin_phi*cos_theta + qc*cos_phi*cos_theta; 
     107 
     108    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
     109    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
     110 
     111    return CALL_IQ_AC(dqa, dqc, table); 
     112} 
     113 
     114#elif defined(CALL_IQ_ABC) // oriented asymmetric 
     115 
     116static double 
     117view_direct(double qx, double qy, 
     118             double theta, double phi, double psi, 
     119             ParameterTable table) 
     120{ 
     121    double sin_theta, cos_theta; 
     122    double sin_phi, cos_phi; 
     123    double sin_psi, cos_psi; 
     124 
     125    // reverse view 
     126    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     127    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     128    SINCOS(psi*M_PI_180, sin_psi, cos_psi); 
     129    const double qa = qx*(sin_phi*sin_psi + cos_phi*cos_psi*cos_theta) + qy*(sin_phi*cos_psi*cos_theta - sin_psi*cos_phi); 
     130    const double qb = qx*(-sin_phi*cos_psi + sin_psi*cos_phi*cos_theta) + qy*(sin_phi*sin_psi*cos_theta + cos_phi*cos_psi); 
     131    const double qc = qx*sin_theta*cos_phi + qy*sin_phi*sin_theta; 
     132 
     133    // reverse jitter after view 
     134    SINCOS(table.theta*M_PI_180, sin_theta, cos_theta); 
     135    SINCOS(table.phi*M_PI_180, sin_phi, cos_phi); 
     136    SINCOS(table.psi*M_PI_180, sin_psi, cos_psi); 
     137    const double dqa = qa*cos_psi*cos_theta + qb*(sin_phi*sin_theta*cos_psi - sin_psi*cos_phi) + qc*(-sin_phi*sin_psi - sin_theta*cos_phi*cos_psi); 
     138    const double dqb = qa*sin_psi*cos_theta + qb*(sin_phi*sin_psi*sin_theta + cos_phi*cos_psi) + qc*(sin_phi*cos_psi - sin_psi*sin_theta*cos_phi); 
     139    const double dqc = qa*sin_theta - qb*sin_phi*cos_theta + qc*cos_phi*cos_theta; 
     140 
     141    return CALL_IQ_ABC(dqa, dqb, dqc, table); 
     142} 
     143#endif 
     144 
     145#endif // !MAGNETIC 
    73146 
    74147kernel 
    … …  
    111184#endif // MAGNETIC 
    112185 
     186#if defined(CALL_IQ_AC) // oriented symmetric 
     187  const double theta = values[details->theta_par+2]; 
     188  const double phi = values[details->theta_par+3]; 
     189#elif defined(CALL_IQ_ABC) // oriented asymmetric 
     190  const double theta = values[details->theta_par+2]; 
     191  const double phi = values[details->theta_par+3]; 
     192  const double psi = values[details->theta_par+4]; 
     193#endif 
     194 
    113195  // Fill in the initial variables 
    114196  //   values[0] is scale 
    … …  
    215297 
    216298//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n"); } 
     299//#if defined(CALL_IQ_AC) // oriented symmetric 
     300//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); printf("theta=%g dtheta=%g",theta,local_values.table.theta); printf("\n"); } 
     301//#endif 
    217302 
    218303    #ifdef INVALID 
    … …  
    267352                } 
    268353                #endif 
    269                 scattering += CALL_IQ(q, q_index, local_values.table); 
     354#  if defined(CALL_IQ_A) // unoriented 
     355                scattering += CALL_IQ_A(sqrt(qsq), local_values.table); 
     356#  elif defined(CALL_IQ_AC) // oriented symmetric 
     357                scattering += view_direct(qx, qy, theta, phi, local_values.table); 
     358#  elif defined(CALL_IQ_ABC) // oriented asymmetric 
     359                scattering += view_direct(qx, qy, theta, phi, psi, local_values.table); 
     360#  endif 
    270361              } 
    271362            } 
    272363          } 
    273364        } 
    274 #else  // !MAGNETIC 
    275         const double scattering = CALL_IQ(q, q_index, local_values.table); 
     365#elif defined(CALL_IQ) // 1d, not MAGNETIC 
     366          const double scattering = CALL_IQ(q[q_index], local_values.table); 
     367#else  // 2d data, not magnetic 
     368          const double qx = q[2*q_index]; 
     369          const double qy = q[2*q_index+1]; 
     370#  if defined(CALL_IQ_A) // unoriented 
     371          const double scattering = CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table); 
     372#  elif defined(CALL_IQ_AC) // oriented symmetric 
     373          const double scattering = view_direct(qx, qy, theta, phi, local_values.table); 
     374#  elif defined(CALL_IQ_ABC) // oriented asymmetric 
     375          const double scattering = view_direct(qx, qy, theta, phi, psi, local_values.table); 
     376#  endif 
    276377#endif // !MAGNETIC 
    277378        this_result += weight0 * scattering; 

  • sasmodels/kernelpy.py

    r62d7601 r8698a0d  
    113113 
    114114        partable = model_info.parameters 
    115         kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 
    116                              else partable.iq_parameters) 
     115        #kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 
     116        #                     else partable.iq_parameters) 
     117        kernel_parameters = partable.iq_parameters 
    117118        volume_parameters = partable.form_volume_parameters 
    118119 

  • sasmodels/modelinfo.py

    ra85a569 r8698a0d  
    382382      with vector parameter p sent as p[]. 
    383383 
    384     * *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 
     384    * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 
    385385      function, with vector parameter p sent as p[]. 
    386386 
    … …  
    420420        # type: (List[Parameter]) -> None 
    421421        self.kernel_parameters = parameters 
     422        self._check_angles() 
    422423        self._set_vector_lengths() 
    423424 
    … …  
    438439        self.iq_parameters = [p for p in self.kernel_parameters 
    439440                              if p.type not in ('orientation', 'magnetic')] 
    440         self.iqxy_parameters = [p for p in self.kernel_parameters 
    441                                 if p.type != 'magnetic'] 
     441        # note: orientation no longer sent to Iqxy, so its the same as 
     442        #self.iqxy_parameters = [p for p in self.kernel_parameters 
     443        #                        if p.type != 'magnetic'] 
    442444        self.form_volume_parameters = [p for p in self.kernel_parameters 
    443445                                       if p.type == 'volume'] 
    … …  
    461463        self.has_2d = any(p.type in ('orientation', 'magnetic') 
    462464                          for p in self.kernel_parameters) 
     465        self.is_asymmetric = any(p.name == 'psi' for p in self.kernel_parameters) 
    463466        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    464467                                if p.id.startswith('M0:')] 
    … …  
    467470                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    468471        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
     472 
     473    def _check_angles(self): 
     474        theta = phi = psi = -1 
     475        for k, p in enumerate(self.kernel_parameters): 
     476            if p.name == 'theta': 
     477                theta = k 
     478                if p.type != 'orientation': 
     479                    raise TypeError("theta must be an orientation parameter") 
     480            elif p.name == 'phi': 
     481                phi = k 
     482                if p.type != 'orientation': 
     483                    raise TypeError("phi must be an orientation parameter") 
     484            elif p.name == 'psi': 
     485                psi = k 
     486                if p.type != 'orientation': 
     487                    raise TypeError("psi must be an orientation parameter") 
     488        if theta >= 0 and phi >= 0: 
     489            if phi != theta+1: 
     490                raise TypeError("phi must follow theta") 
     491            if psi >= 0 and psi != phi+1: 
     492                raise TypeError("psi must follow phi") 
     493        elif theta >= 0 or phi >= 0 or psi >= 0: 
     494            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
    469495 
    470496    def __getitem__(self, key): 

  • sasmodels/weights.py

    rd4c33d6 r8698a0d  
    5757        if not relative: 
    5858            # For orientation, the jitter is relative to 0 not the angle 
    59             #center = 0 
     59            center = 0 
    6060            pass 
    6161        if sigma == 0 or self.npts < 2: 
    … …  
    229229    obj = cls(n, width, nsigmas) 
    230230    v, w = obj.get_weights(value, limits[0], limits[1], relative) 
    231     return v, w 
     231    return v, w/np.sum(w) 
    232232 
    233233 
Note: See TracChangeset for help on using the changeset viewer.