Changeset 32e3c9b in sasmodels


Ignore:
Timestamp:
Jul 21, 2016 4:08:04 PM (8 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
b966a96
Parents:
42356c8
Message:

dll version of magnetic sld

Files:
1 added
1 deleted
13 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/__init__.py

    r59d94b3 r32e3c9b  
    4141    # can build model docs on the fly, including images. 
    4242    return_list = [ 
    43         expand_patterns([], ['*.c', '*.cl', 'convert.json']), 
     43        expand_patterns([], ['*.c', '*.cl']), 
    4444        expand_patterns(['models'], ['*.py', '*.c']), 
    4545        expand_patterns(['models', 'lib'], ['*.c']), 

  • sasmodels/convert.py

    ra0494e9 r32e3c9b  
    77import math 
    88import warnings 
    9 import json 
     9 
     10from .conversion_table import CONVERSION_TABLE 
    1011 
    1112# List of models which SasView versions don't contain the explicit 'scale' argument. 
    … …  
    4748    ] 
    4849 
    49 CONVERSION_TABLE = None 
    50  
    51 def _read_conversion_table(): 
    52     global CONVERSION_TABLE 
    53     if CONVERSION_TABLE is None: 
    54         path = joinpath(dirname(abspath(__file__)), "convert.json") 
    55         with open(path) as fid: 
    56             CONVERSION_TABLE = json_load_byteified(fid) 
    57  
    58 def json_load_byteified(file_handle): 
    59     return _byteify( 
    60         json.load(file_handle, object_hook=_byteify), 
    61         ignore_dicts=True 
    62     ) 
    63  
    64 def _byteify(data, ignore_dicts = False): 
    65     # if this is a unicode string, return its string representation 
    66     if isinstance(data, unicode): 
    67         return data.encode('utf-8') 
    68     # if this is a list of values, return list of byteified values 
    69     if isinstance(data, list): 
    70         return [ _byteify(item, ignore_dicts=True) for item in data ] 
    71     # if this is a dictionary, return dictionary of byteified keys and values 
    72     # but only if we haven't already byteified it 
    73     if isinstance(data, dict) and not ignore_dicts: 
    74         return dict((_byteify(key, ignore_dicts=True), 
    75                      _byteify(value, ignore_dicts=True)) 
    76                     for key, value in data.items()) 
    77     # if it's anything else, return it in its original form 
    78     return data 
    79  
    80  
    8150def _convert_pars(pars, mapping): 
    8251    """ 
    … …  
    9362    return newpars 
    9463 
    95 def _rescale_sld(pars): 
    96     """ 
    97     rescale all sld parameters in the new model definition by 1e6 so the 
    98     numbers are nicer.  Relies on the fact that all sld parameters in the 
    99     new model definition end with sld. 
    100     """ 
    101     return dict((p, (v*1e6 if p.endswith('sld') 
    102                      else v*1e-15 if 'ndensity' in p 
    103                      else v)) 
    104                 for p, v in pars.items()) 
    105  
    10664def convert_model(name, pars): 
    10765    """ 
    … …  
    11674    return [pk*scale for pk in par] if isinstance(par, list) else par*scale 
    11775 
    118 def _unscale_sld(pars): 
     76def _is_sld(modelinfo, id): 
     77    if id.startswith('M0:'): 
     78        return True 
     79    if (id.endswith('_pd') or id.endswith('_pd_n') or id.endswith('_pd_nsigma') 
     80            or id.endswith('_pd_width') or id.endswith('_pd_type')): 
     81        return False 
     82    for p in modelinfo.parameters.call_parameters: 
     83        if p.id == id: 
     84            return p.type == 'sld' 
     85    # check through kernel parameters in case it is a named as a vector 
     86    for p in modelinfo.parameters.kernel_parameters: 
     87        if p.id == id: 
     88            return p.type == 'sld' 
     89    raise ValueError("unknown parameter %r in conversion"%id) 
     90 
     91def _unscale_sld(modelinfo, pars): 
    11992    """ 
    12093    rescale all sld parameters in the new model definition by 1e6 so the 
    … …  
    12295    new model definition end with sld. 
    12396    """ 
    124     return dict((p, (_unscale(v,1e-6) if p.startswith('sld') or p.endswith('sld') 
    125                      else _unscale(v,1e15) if 'ndensity' in p 
    126                      else v)) 
    127                 for p, v in pars.items()) 
     97    return dict((id, (_unscale(v,1e-6) if _is_sld(modelinfo, id) else v)) 
     98                for id, v in pars.items()) 
    12899 
    129100def _remove_pd(pars, key, name): 
    … …  
    164135 
    165136def revert_name(model_info): 
    166     _read_conversion_table() 
    167137    oldname, oldpars = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    168138    return oldname 
    169139 
    170140def _get_translation_table(model_info): 
    171     _read_conversion_table() 
    172141    _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    173142    translation = translation.copy() 
    … …  
    186155 
    187156def _trim_vectors(model_info, pars, oldpars): 
    188     _read_conversion_table() 
    189157    _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    190158    for p in model_info.parameters.kernel_parameters: 
    … …  
    212180    else: 
    213181        translation = _get_translation_table(model_info) 
    214         oldpars = _revert_pars(_unscale_sld(pars), translation) 
     182        oldpars = _revert_pars(_unscale_sld(model_info, pars), translation) 
    215183        oldpars = _trim_vectors(model_info, pars, oldpars) 
    216184 
    … …  
    243211            _remove_pd(oldpars, 'rimB', name) 
    244212            _remove_pd(oldpars, 'rimC', name) 
     213        elif name == 'polymer_micelle': 
     214            if 'ndensity' in oldpars: 
     215                oldpars['ndensity'] *= 1e15 
    245216        elif name == 'spherical_sld': 
    246217            for k in range(1, int(pars['n_shells'])+1): 

  • sasmodels/details.py

    ra738209 r32e3c9b  
    9696 
    9797    pd_offset = np.cumsum(np.hstack((0, pd_length))) 
     98    #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(model_info.parameters.call_parameters))) 
     99    #print("len:",pd_length) 
     100    #print("off:",pd_offset) 
    98101    # Note: the reversing view, x[::-1], does not require a copy 
    99102    idx = np.argsort(pd_length)[::-1][:num_active] 

  • sasmodels/direct_model.py

    ra738209 r32e3c9b  
    7070 
    7171    call_details, values = kernel.build_details(calculator, vw_pairs) 
    72     return calculator(call_details, values, cutoff) 
     72    magnetic = any(values[k]!=0 for k in parameters.magnetism_index) 
     73    #print("values:", values) 
     74    return calculator(call_details, values, cutoff, magnetic) 
    7375 
    7476def get_weights(parameter, values): 

  • sasmodels/generate.py

    ra738209 r32e3c9b  
    387387 
    388388 
    389 def kernel_name(model_info, is_2d): 
    390     # type: (ModelInfo, bool) -> str 
     389def kernel_name(model_info, variant): 
     390    # type: (ModelInfo, str) -> str 
    391391    """ 
    392392    Name of the exported kernel symbol. 
    393     """ 
    394     return model_info.name + "_" + ("Iqxy" if is_2d else "Iq") 
     393 
     394    *variant* is "Iq", "Iqxy" or "Imagnetic". 
     395    """ 
     396    return model_info.name + "_" + variant 
    395397 
    396398 
    … …  
    555557    refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 
    556558    call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 
    557     if _have_Iqxy(user_code): 
     559    if _have_Iqxy(user_code) or isinstance(model_info.Iqxy, str): 
    558560        # Call 2D model 
    559         refs = ["q[2*_i]", "q[2*_i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 
     561        refs = ["_q[2*_i]", "_q[2*_i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 
    560562        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 
    561563    else: 
    … …  
    566568        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
    567569 
     570    magpars = [k-2 for k,p in enumerate(partable.call_parameters) 
     571               if p.type == 'sld'] 
     572 
    568573    # Fill in definitions for numbers of parameters 
    569574    source.append("#define MAX_PD %s"%partable.max_pd) 
    570575    source.append("#define NPARS %d"%partable.npars) 
     576    source.append("#define NUM_MAGNETIC %d" % len(magpars)) 
     577    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    571578 
    572579    # TODO: allow mixed python/opencl kernels? 
    … …  
    584591    source = [ 
    585592        # define the Iq kernel 
    586         "#define KERNEL_NAME %s_Iq"%name, 
     593        "#define KERNEL_NAME %s_Iq" % name, 
    587594        call_iq, 
    588595        kernel_code, 
    … …  
    591598 
    592599        # define the Iqxy kernel from the same source with different #defines 
    593         "#define KERNEL_NAME %s_Iqxy"%name, 
     600        "#define KERNEL_NAME %s_Iqxy" % name, 
    594601        call_iqxy, 
    595602        kernel_code, 
     603        "#undef CALL_IQ", 
     604        "#undef KERNEL_NAME", 
     605 
     606        # define the Imagnetic kernel 
     607        "#define KERNEL_NAME %s_Imagnetic" % name, 
     608        "#define MAGNETIC 1", 
     609        call_iqxy, 
     610        kernel_code, 
     611        "#undef MAGNETIC", 
    596612        "#undef CALL_IQ", 
    597613        "#undef KERNEL_NAME", 

  • sasmodels/kernel.py

    rdef2c1b r32e3c9b  
    4141    results = None # type: List[np.ndarray] 
    4242 
    43     def __call__(self, call_details, values, cutoff): 
    44         # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 
     43    def __call__(self, call_details, values, cutoff, magnetic): 
     44        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    4545        raise NotImplementedError("need to implement __call__") 
    4646 

  • sasmodels/kernel_iq.c

    ra738209 r32e3c9b  
    3333#endif 
    3434 
     35#ifdef MAGNETIC 
     36const int32_t magnetic[] = { MAGNETIC_PARS }; 
     37#endif 
     38 
     39#ifdef MAGNETIC 
     40// Return value restricted between low and high 
     41static double clip(double value, double low, double high) 
     42{ 
     43    return (value < low ? low : (value > high ? high : value)); 
     44} 
     45 
     46// Compute spin cross sections given in_spin and out_spin 
     47// To convert spin cross sections to sld b: 
     48//     uu * (sld - m_sigma_x); 
     49//     dd * (sld + m_sigma_x); 
     50//     ud * (m_sigma_y + 1j*m_sigma_z); 
     51//     du * (m_sigma_y - 1j*m_sigma_z); 
     52static void spins(double in_spin, double out_spin, 
     53    double *uu, double *dd, double *ud, double *du) 
     54{ 
     55    in_spin = clip(in_spin, 0.0, 1.0); 
     56    out_spin = clip(out_spin, 0.0, 1.0); 
     57        *uu = sqrt(sqrt(in_spin * out_spin)); 
     58        *dd = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); 
     59        *ud = sqrt(sqrt(in_spin * (1.0-out_spin))); 
     60        *du = sqrt(sqrt((1.0-in_spin) * out_spin)); 
     61} 
     62 
     63// Convert polar to rectangular coordinates. 
     64static void polrec(double r, double theta, double phi, 
     65    double *x, double *y, double *z) 
     66{ 
     67    double cos_theta, sin_theta, cos_phi, sin_phi; 
     68    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     69    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     70    *x = r * cos_theta * cos_phi; 
     71    *y = r * sin_theta; 
     72    *z = -r * cos_theta * sin_phi; 
     73} 
     74#endif 
    3575 
    3676kernel 
    … …  
    4080    const int32_t pd_stop,      // where we are stopping in the polydispersity loop 
    4181    global const ProblemDetails *details, 
    42     global const double *values, 
     82   // global const  // TODO: make it const again! 
     83    double *values, 
    4384    global const double *q, // nq q values, with padding to boundary 
    4485    global double *result,  // nq+3 return values, again with padding 
    … …  
    5293 
    5394  // Fill in the initial variables 
     95  //   values[0] is scale 
     96  //   values[1] is background 
    5497  #ifdef USE_OPENMP 
    5598  #pragma omp parallel for 
    … …  
    58101    pvec[k] = values[k+2]; 
    59102  } 
     103#ifdef MAGNETIC 
     104  const double up_frac_i = values[NPARS+2]; 
     105  const double up_frac_f = values[NPARS+3]; 
     106  const double up_angle = values[NPARS+4]; 
     107  #define MX(_k) (values[NPARS+5+3*_k]) 
     108  #define MY(_k) (values[NPARS+6+3*_k]) 
     109  #define MZ(_k) (values[NPARS+7+3*_k]) 
     110 
     111  // TODO: precompute this on the python side 
     112  // Convert polar to rectangular coordinates in place. 
     113  if (pd_start == 0) {  // Update in place; only do this for the first hunk! 
     114//printf("spin: %g %g %g\n", up_frac_i, up_frac_f, up_angle); 
     115    for (int mag=0; mag < NUM_MAGNETIC; mag++) { 
     116//printf("mag %d: %g %g %g\n", mag, MX(mag), MY(mag), MZ(mag)); 
     117        polrec(MX(mag), MY(mag), MZ(mag), &MX(mag), &MY(mag), &MZ(mag)); 
     118//printf("   ==>: %g %g %g\n", MX(mag), MY(mag), MZ(mag)); 
     119    } 
     120  } 
     121  // Interpret polarization cross section. 
     122  double uu, dd, ud, du; 
     123  double cos_mspin, sin_mspin; 
     124  spins(up_frac_i, up_frac_f, &uu, &dd, &ud, &du); 
     125  SINCOS(-up_angle*M_PI_180, sin_mspin, cos_mspin); 
     126#endif 
    60127 
    61128  // Monodisperse computation 
    … …  
    74141    #endif 
    75142    for (int q_index=0; q_index < nq; q_index++) { 
     143#ifdef MAGNETIC 
     144      const double qx = q[2*q_index]; 
     145      const double qy = q[2*q_index+1]; 
     146      const double qsq = qx*qx + qy*qy; 
     147 
     148      // Constant across orientation, polydispersity for given qx, qy 
     149      double px, py, pz; 
     150      if (qsq > 1e-16) { 
     151        px = (qy*cos_mspin + qx*sin_mspin)/qsq; 
     152        py = (qy*sin_mspin - qx*cos_mspin)/qsq; 
     153        pz = 1.0; 
     154      } else { 
     155        px = py = pz = 0.0; 
     156      } 
     157 
     158      double scattering = 0.0; 
     159      if (uu > 1e-8) { 
     160        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     161            const double perp = (qy*MX(mag) - qx*MY(mag)); 
     162            pvec[magnetic[mag]] = (values[magnetic[mag]+2] - perp*px)*uu; 
     163        } 
     164        scattering += CALL_IQ(q, q_index, local_values); 
     165      } 
     166      if (dd > 1e-8){ 
     167        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     168            const double perp = (qy*MX(mag) - qx*MY(mag)); 
     169            pvec[magnetic[mag]] = (values[magnetic[mag]+2] + perp*px)*dd; 
     170        } 
     171        scattering += CALL_IQ(q, q_index, local_values); 
     172      } 
     173      if (ud > 1e-8){ 
     174        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     175            const double perp = (qy*MX(mag) - qx*MY(mag)); 
     176            pvec[magnetic[mag]] = perp*py*ud; 
     177        } 
     178        scattering += CALL_IQ(q, q_index, local_values); 
     179        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     180            pvec[magnetic[mag]] = MZ(mag)*pz*ud; 
     181        } 
     182        scattering += CALL_IQ(q, q_index, local_values); 
     183      } 
     184      if (du > 1e-8) { 
     185        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     186            const double perp = (qy*MX(mag) - qx*MY(mag)); 
     187            pvec[magnetic[mag]] = perp*py*du; 
     188        } 
     189        scattering += CALL_IQ(q, q_index, local_values); 
     190        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     191            pvec[magnetic[mag]] = -MZ(mag)*pz*du; 
     192        } 
     193        scattering += CALL_IQ(q, q_index, local_values); 
     194      } 
     195#else 
    76196      double scattering = CALL_IQ(q, q_index, local_values); 
     197#endif 
    77198      result[q_index] = (norm>0. ? scale*scattering/norm + background : background); 
    78199    } 
    … …  
    82203#if MAX_PD > 0 
    83204 
     205#if MAGNETIC 
     206  const double *pd_value = values+2+NPARS+3+3*NUM_MAGNETIC; 
     207#else 
    84208  const double *pd_value = values+2+NPARS; 
     209#endif 
    85210  const double *pd_weight = pd_value+details->pd_sum; 
    86211 
    … …  
    164289      #endif 
    165290      for (int q_index=0; q_index < nq; q_index++) { 
    166         const double scattering = CALL_IQ(q, q_index, local_values); 
     291#ifdef MAGNETIC 
     292        const double qx = q[2*q_index]; 
     293        const double qy = q[2*q_index+1]; 
     294        const double qsq = qx*qx + qy*qy; 
     295 
     296        // Constant across orientation, polydispersity for given qx, qy 
     297        double px, py, pz; 
     298        if (qsq > 1e-16) { 
     299          px = (qy*cos_mspin + qx*sin_mspin)/qsq; 
     300          py = (qy*sin_mspin - qx*cos_mspin)/qsq; 
     301          pz = 1.0; 
     302        } else { 
     303          px = py = pz = 0.0; 
     304        } 
     305 
     306        double scattering = 0.0; 
     307        if (uu > 1e-8) { 
     308          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     309              const double perp = (qy*MX(mag) - qx*MY(mag)); 
     310              pvec[magnetic[mag]] = (values[magnetic[mag]+2] - perp*px)*uu; 
     311          } 
     312          scattering += CALL_IQ(q, q_index, local_values); 
     313        } 
     314        if (dd > 1e-8){ 
     315          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     316              const double perp = (qy*MX(mag) - qx*MY(mag)); 
     317              pvec[magnetic[mag]] = (values[magnetic[mag]+2] + perp*px)*dd; 
     318          } 
     319          scattering += CALL_IQ(q, q_index, local_values); 
     320        } 
     321        if (ud > 1e-8){ 
     322          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     323              const double perp = (qy*MX(mag) - qx*MY(mag)); 
     324              pvec[magnetic[mag]] = perp*py*ud; 
     325          } 
     326          scattering += CALL_IQ(q, q_index, local_values); 
     327          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     328              pvec[magnetic[mag]] = MZ(mag)*pz*ud; 
     329          } 
     330          scattering += CALL_IQ(q, q_index, local_values); 
     331        } 
     332        if (du > 1e-8) { 
     333          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     334              const double perp = (qy*MX(mag) - qx*MY(mag)); 
     335              pvec[magnetic[mag]] = perp*py*du; 
     336          } 
     337          scattering += CALL_IQ(q, q_index, local_values); 
     338          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     339              pvec[magnetic[mag]] = -MZ(mag)*pz*du; 
     340          } 
     341          scattering += CALL_IQ(q, q_index, local_values); 
     342        } 
     343#else 
     344        double scattering = CALL_IQ(q, q_index, local_values); 
     345#endif 
    167346        result[q_index] += weight*scattering; 
    168347      } 

  • sasmodels/kernelcl.py

    rdef2c1b r32e3c9b  
    385385        self.program = None 
    386386 
    387     def make_kernel(self, q_vectors): 
     387    def make_kernel(self, q_vectors, magnetic=False): 
    388388        # type: (List[np.ndarray]) -> "GpuKernel" 
    389389        if self.program is None: 
    … …  
    392392                                    self.dtype, self.fast) 
    393393        is_2d = len(q_vectors) == 2 
    394         kernel_name = generate.kernel_name(self.info, is_2d) 
     394        variant = "Imagnetic" if magnetic else "Iqxy" if is_2d else "Iq" 
     395        kernel_name = generate.kernel_name(self.info, variant) 
    395396        kernel = getattr(self.program, kernel_name) 
    396397        return GpuKernel(kernel, self.dtype, self.info, q_vectors) 
    … …  
    519520                     else np.float32)  # will never get here, so use np.float32 
    520521 
    521     def __call__(self, call_details, values, cutoff): 
    522         # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 
     522    def __call__(self, call_details, values, cutoff, magnetic): 
     523        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    523524        context = self.queue.context 
    524525        # Arrange data transfer to card 

  • sasmodels/kerneldll.py

    rdef2c1b r32e3c9b  
    270270        # int, int, int, int*, double*, double*, double*, double*, double 
    271271        argtypes = [c_int32]*3 + [c_void_p]*4 + [fp] 
    272         self._Iq = self._dll[generate.kernel_name(self.info, is_2d=False)] 
    273         self._Iqxy = self._dll[generate.kernel_name(self.info, is_2d=True)] 
     272        self._Iq = self._dll[generate.kernel_name(self.info, "Iq")] 
     273        self._Iqxy = self._dll[generate.kernel_name(self.info, "Iqxy")] 
     274        self._Imagnetic = self._dll[generate.kernel_name(self.info, "Imagnetic")] 
    274275        self._Iq.argtypes = argtypes 
    275276        self._Iqxy.argtypes = argtypes 
     277        self._Imagnetic.argtypes = argtypes 
    276278 
    277279    def __getstate__(self): 
    … …  
    284286        self._dll = None 
    285287 
    286     def make_kernel(self, q_vectors): 
     288    def make_kernel(self, q_vectors, magnetic=False): 
    287289        # type: (List[np.ndarray]) -> DllKernel 
    288290        q_input = PyInput(q_vectors, self.dtype) 
    … …  
    290292        if self._dll is None: 
    291293            self._load_dll() 
    292         kernel = self._Iqxy if q_input.is_2d else self._Iq 
     294        kernel = [self._Iqxy, self._Imagnetic] if q_input.is_2d else self._Iq 
    293295        return DllKernel(kernel, self.info, q_input) 
    294296 
    … …  
    341343                     else np.float128) 
    342344 
    343     def __call__(self, call_details, values, cutoff): 
    344         # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 
     345    def __call__(self, call_details, values, cutoff, magnetic): 
     346        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    345347 
    346348        #print("in kerneldll") 
    … …  
    358360            ] 
    359361        #print("calling DLL") 
    360         self.kernel(*args) # type: ignore 
     362        self.kernel[1 if magnetic else 0](*args) # type: ignore 
    361363        return self.result[:-1] 
    362364 

  • sasmodels/kernelpy.py

    ra738209 r32e3c9b  
    155155                        else (lambda: 1.0)) 
    156156 
    157     def __call__(self, call_details, values, cutoff): 
     157    def __call__(self, call_details, values, cutoff, magnetic): 
    158158        assert isinstance(call_details, details.CallDetails) 
    159159        res = _loops(self._parameter_vector, self._form, self._volume, 

  • sasmodels/modelinfo.py

    rce176ca r32e3c9b  
    414414        self.kernel_parameters = parameters 
    415415        self._set_vector_lengths() 
     416 
     417        self.npars = sum(p.length for p in self.kernel_parameters) 
     418        self.num_magnetic = sum(p.length for p in self.kernel_parameters 
     419                                if p.type=='sld') 
     420 
    416421        self.call_parameters = self._get_call_parameters() 
    417422        self.defaults = self._get_defaults() 
    … …  
    441446        num_pd = sum(p.length for p in self.kernel_parameters if p.polydisperse) 
    442447        # Don't use more polydisperse parameters than are available in the model 
    443         # Note: we can do polydispersity on arbitrary parameters, so it is not 
    444         # clear that this is a good idea; it does however make the poly_details 
    445         # code easier to write, so we will leave it in for now. 
    446448        self.max_pd = min(num_pd, MAX_PD) 
    447  
    448         self.npars = sum(p.length for p in self.kernel_parameters) 
    449449 
    450450        # true if has 2D parameters 
    451451        self.has_2d = any(p.type in ('orientation', 'magnetic') 
    452452                          for p in self.kernel_parameters) 
     453        self.magnetism_index = [k for k,p in enumerate(self.call_parameters) 
     454                                if p.id.startswith('M0:')] 
    453455 
    454456        self.pd_1d = set(p.name for p in self.call_parameters 
    455457                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    456         self.pd_2d = set(p.name for p in self.call_parameters 
    457                          if p.polydisperse and p.type != 'magnetic') 
     458        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
    458459 
    459460    def _set_vector_lengths(self): 
    … …  
    519520                    pk.relative_pd = p.relative_pd 
    520521                    full_list.append(pk) 
     522 
     523        # Add the magnetic parameters to the end of the call parameter list. 
     524        if self.num_magnetic > 0: 
     525            full_list.extend([ 
     526                Parameter('up:frac_i', '', 0., [0., 1.], 
     527                          'magnetic', 'fraction of spin up incident'), 
     528                Parameter('up:frac_f', '', 0., [0., 1.], 
     529                          'magnetic', 'fraction of spin up final'), 
     530                Parameter('up:angle', 'degress', 0., [0., 360.], 
     531                          'magnetic', 'spin up angle'), 
     532            ]) 
     533            slds = [p for p in full_list if p.type == 'sld'] 
     534            for p in slds: 
     535                full_list.extend([ 
     536                    Parameter('M0:'+p.id, '1e-6/Ang^2', 0., [-np.inf, np.inf], 
     537                              'magnetic', 'magnetic amplitude for '+p.description), 
     538                    Parameter('mtheta:'+p.id, 'degrees', 0., [-90., 90.], 
     539                               'magnetic', 'magnetic latitude for '+p.description), 
     540                    Parameter('mphi:'+p.id, 'degrees', 0., [-180., 180.], 
     541                               'magnetic', 'magnetic longitude for '+p.description), 
     542                ]) 
     543        #print("call parameters", full_list) 
    521544        return full_list 
    522545 
    … …  
    553576        parameter is always returned first since the GUI will want to set it 
    554577        early, and rerender the table when it is changed. 
     578 
     579        Parameters marked as sld will automatically have a set of associated 
     580        magnetic parameters (m0:p, mtheta:p, mphi:p), as well as polarization 
     581        information (up:theta, up:frac_i, up:frac_f). 
    555582        """ 
    556583        # control parameters go first 
    … …  
    565592 
    566593        # Gather entries such as name[4] into groups of the same length 
    567         fixed = {}  # type: Dict[int, List[Parameter]] 
     594        fixed_length = {}  # type: Dict[int, List[Parameter]] 
    568595        for p in self.kernel_parameters: 
    569596            if p.length > 1 and p.length_control is None: 
    570                 fixed.setdefault(p.length, []).append(p) 
     597                fixed_length.setdefault(p.length, []).append(p) 
    571598 
    572599        # Using the call_parameters table, we already have expanded forms 
    573600        # for each of the vector parameters; put them in a lookup table 
    574601        expanded_pars = dict((p.name, p) for p in self.call_parameters) 
     602 
     603        def append_group(name): 
     604            """add the named parameter, and related magnetic parameters if any""" 
     605            result.append(expanded_pars[name]) 
     606            if is2d: 
     607                for tag in 'M0:', 'mtheta:', 'mphi:': 
     608                    if tag+name in expanded_pars: 
     609                        result.append(expanded_pars[tag+name]) 
    575610 
    576611        # Gather the user parameters in order 
    … …  
    589624                    for k in range(1, table_length+1): 
    590625                        for entry in table: 
    591                             result.append(expanded_pars[entry.id+str(k)]) 
     626                            append_group(entry.id+str(k)) 
    592627                else: 
    593628                    pass # already processed all entries 
    594629            elif p.length > 1: 
    595                 table = fixed.get(p.length, []) 
     630                table = fixed_length.get(p.length, []) 
    596631                if table: 
    597632                    table_length = p.length 
    598                     del fixed[p.length] 
     633                    del fixed_length[p.length] 
    599634                    for k in range(1, table_length+1): 
    600635                        for entry in table: 
    601                             result.append(expanded_pars[entry.id+str(k)]) 
     636                            append_group(entry.id+str(k)) 
    602637                else: 
    603638                    pass # already processed all entries 
    604639            else: 
    605                 result.append(p) 
     640                append_group(p.id) 
     641 
     642        if is2d and 'up:angle' in expanded_pars: 
     643            result.extend([ 
     644                expanded_pars['up:frac_i'], 
     645                expanded_pars['up:frac_f'], 
     646                expanded_pars['up:angle'], 
     647            ]) 
    606648 
    607649        return result 
    … …  
    634676    if (model_info.Iq is None 
    635677        and model_info.Iqxy is None 
     678        and model_info.Imagnetic is None 
    636679        and model_info.form_volume is None): 
    637680        return 
    … …  
    643686        return 
    644687    for k, v in enumerate(source.split('\n')): 
    645         if v.startswith('Iqxy'): 
     688        if v.startswith('Imagnetic'): 
     689            model_info._Imagnetic_line = k+1 
     690        elif v.startswith('Iqxy'): 
    646691            model_info._Iqxy_line = k+1 
    647692        elif v.startswith('Iq'): 
    … …  
    692737    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    693738    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 
     739    info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore 
    694740    info.profile = getattr(kernel_module, 'profile', None) # type: ignore 
    695741    info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore 
    … …  
    843889    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
    844890    Iqxy = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     891    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
     892    Imagnetic = None        # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    845893    #: Returns a model profile curve *x, y*.  If *profile* is defined, this 
    846894    #: curve will appear in response to the *Show* button in SasView.  Use 
    … …  
    857905    # line numbers within the python file for bits of C source, if defined 
    858906    # NB: some compilers fail with a "#line 0" directive, so default to 1. 
     907    _Imagnetic_line = 1 
    859908    _Iqxy_line = 1 
    860909    _Iq_line = 1 

  • sasmodels/sasview_model.py

    rce176ca r32e3c9b  
    152152            orientation_params.append(p.name+".width") 
    153153            fixed.append(p.name+".width") 
    154         if p.type == 'magnetic': 
     154        elif p.type == 'magnetic': 
    155155            orientation_params.append(p.name) 
    156156            magnetic_params.append(p.name) 
    157157            fixed.append(p.name+".width") 
     158 
    158159 
    159160    # Build class dictionary 
    … …  
    512513                 for p in self._model_info.parameters.call_parameters] 
    513514        call_details, value = kernel.build_details(calculator, pairs) 
    514         result = calculator(call_details, value, cutoff=self.cutoff) 
     515        result = calculator(call_details, value, cutoff=self.cutoff, magnetic=False) 
    515516        calculator.release() 
    516517        return result 

  • setup.py

    r3a161d0 r32e3c9b  
    4646    package_data={ 
    4747        'sasmodels.models': ['*.c', 'lib/*.c'], 
    48         'sasmodels': ['*.c', '*.cl', 'convert.json'], 
     48        'sasmodels': ['*.c', '*.cl'], 
    4949    }, 
    5050    install_requires = [ 
Note: See TracChangeset for help on using the changeset viewer.