Changeset a7c2cfe in sasmodels


Ignore:
Timestamp:
Apr 18, 2016 9:47:24 AM (9 years ago)
Author:
wojciech
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:
1bf66d9
Parents:
4605bf10 (diff), 38a9b07 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'polydisp' of ​https://github.com/SasView/sasmodels into polydisp

Location:
sasmodels
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/model_test.py

    r0ff62d4 r38a9b07  
    259259                    self.assertTrue(np.isfinite(actual_yi), 
    260260                                    'invalid f(%s): %s' % (xi, actual_yi)) 
     261                elif np.isnan(yi): 
     262                    self.assertTrue(np.isnan(actual_yi), 
     263                                    'f(%s): expected:%s; actual:%s' 
     264                                    % (xi, yi, actual_yi)) 
    261265                else: 
    262                     self.assertTrue(is_near(yi, actual_yi, 5), 
     266                    # is_near does not work for infinite values, so also test 
     267                    # for exact values.  Note that this will not 
     268                    self.assertTrue(yi==actual_yi or is_near(yi, actual_yi, 5), 
    263269                                    'f(%s); expected:%s; actual:%s' 
    264270                                    % (xi, yi, actual_yi)) 

  • sasmodels/models/__init__.py

    r32c160a r1ca1fd9  
     1""" 
     2    1D Modeling for SAS 
     3""" 
     4#from sas.models import * 
     5 
     6import os 
     7from distutils.filelist import findall 
     8 
     9__version__ = "2.1.0" 
     10 
     11def get_data_path(media): 
     12    """ 
     13    """ 
     14    # Check for data path in the package 
     15    path = os.path.join(os.path.dirname(__file__), media) 
     16    if os.path.isdir(path): 
     17        return path 
     18 
     19    # Check for data path next to exe/zip file. 
     20    # If we are inside a py2exe zip file, we need to go up 
     21    # to get to the directory containing  
     22    # the media for this module 
     23    path = os.path.dirname(__file__) 
     24    #Look for maximum n_dir up of the current dir to find media 
     25    n_dir = 12 
     26    for i in range(n_dir): 
     27        path, _ = os.path.split(path) 
     28        media_path = os.path.join(path, media) 
     29        if os.path.isdir(media_path): 
     30            module_media_path = os.path.join(media_path, 'models_media') 
     31            if os.path.isdir(module_media_path): 
     32                return module_media_path 
     33            return media_path 
     34    
     35    raise RuntimeError('Could not find models media files') 
     36 
     37def data_files(): 
     38    """ 
     39    Return the data files associated with media. 
     40     
     41    The format is a list of (directory, [files...]) pairs which can be 
     42    used directly in setup(...,data_files=...) for setup.py. 
     43 
     44    """ 
     45    data_files = [] 
     46    path_img = get_data_path(media=os.path.join("sasmodels","sasmodels","models","img")) 
     47    #path_img = get_data_path(media="img") 
     48    im_list = findall(path_img) 
     49    #for f in findall(path): 
     50    #    if os.path.isfile(f) and f not in im_list: 
     51    #        data_files.append(('media/models_media', [f])) 
     52     
     53    for f in im_list: 
     54        data_files.append(('media/models_media/img', [f])) 
     55    return data_files 
     56 

  • sasmodels/models/porod.py

    rec45c4f r82923a6  
    2626""" 
    2727 
    28 from numpy import sqrt, power 
     28from numpy import sqrt, power, inf, errstate 
    2929 
    3030name = "porod" 
    … …  
    4242    @param q: Input q-value 
    4343    """ 
    44     return 1.0/power(q, 4) 
     44    with errstate(divide='ignore'): 
     45        return power(q, -4) 
    4546 
    4647Iq.vectorized = True  # Iq accepts an array of q values 
    … …  
    5859demo = dict(scale=1.5, background=0.5) 
    5960 
    60 tests = [[{'scale': 0.00001, 'background':0.01}, 0.04, 3.916250]] 
     61tests = [ 
     62    [{'scale': 0.00001, 'background':0.01}, 0.04, 3.916250], 
     63    [{}, 0.0, inf], 
     64] 

  • sasmodels/models/spherical_sld.c

    r299dcce r4605bf10  
    1 //Headers 
    2 static double form_volume(double thick_inter[], 
    3     double thick_flat_[], 
    4     double core_radius, 
    5     int n_shells); 
    6  
    7 double Iq(double q, 
     1static double form_volume( 
    82    int n_shells, 
    9     double thick_inter[], 
    10     double func_inter[], 
    11     double sld_core, 
    12     double sld_solvent, 
    13     double sld_flat[], 
     3    double radius_core, 
    144    double thick_flat[], 
    15     double nu_inter[], 
    16     int npts_inter, 
    17     double core_radius); 
    18  
    19 double Iqxy(double qx, double qy, 
    20     int n_shells, 
    21     double thick_inter[], 
    22     double func_inter[], 
    23     double sld_core, 
    24     double sld_solvent, 
    25     double sld_flat[], 
    26     double thick_flat[], 
    27     double nu_inter[], 
    28     int npts_inter, 
    29     double core_radius); 
    30  
    31 //Main code 
    32 static double form_volume(double thick_inter[], 
    33     double thick_flat_[], 
    34     double core_radius, 
    35     int n) 
     5    double thick_inter[]) 
    366{ 
    377    double radius = 0.0; 
    388    int i; 
    39     double r = core_radius; 
    40     for (i=0; i < n; i++) { 
     9    double r = radius_core; 
     10    for (i=0; i < n_shells; i++) { 
    4111        r += thick_inter[i]; 
    4212        r += thick_flat[i]; 
    … …  
    5626  double background = dp[6]; 
    5727  double npts = dp[57]; //number of sub_layers in each interface 
    58   double nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now 
     28  double nsl=npts;//21.0; //nsl = Num_sub_layer:  must be ODD double number 
    5929  int n_s; 
    6030 
    … …  
    6333  double pi; 
    6434 
    65   //int* fun_type; 
    66   //double* sld; 
    67   //double* thick_inter; 
    68   //double* thick; 
    69   //double* fun_coef; 
    70  
    7135  double total_thick=0.0; 
    7236 
    73   //fun_type = (int*)malloc((n+2)*sizeof(int)); 
    74   //sld = (double*)malloc((n+2)*sizeof(double)); 
    75   //thick_inter = (double*)malloc((n+2)*sizeof(double)); 
    76   //thick = (double*)malloc((n+2)*sizeof(double)); 
    77   //fun_coef = (double*)malloc((n+2)*sizeof(double)); 
    78  
    79   //TODO: Solution to avoid mallocs but probablyu can be done better 
    8037  int fun_type[12]; 
    8138  double sld[12]; 
    … …  
    175132            if (fabs(slope) > 0.0 ){ 
    176133              //fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr); 
    177               fun = sign * r * sph_j1c(qr)  +  sign * 3.0 * sin(qr)/(qr * qr * q ) + sign * 6.0 * cos(qr)/(qr * qr * qr * q); 
     134              fun = sign * r * sph_j1c(qr) + sign * 3.0 * sin(qr)/(qr * qr * q ) 
     135                + sign * 6.0 * cos(qr)/(qr * qr * qr * q); 
    178136            } 
    179137          } 
    … …  
    203161    } 
    204162  } 
    205   //vol += vol_sub; 
    206163  f2 = f * f / vol; 
    207   //f2 *= scale; 
    208   //f2 += background; 
    209   //free(fun_type); 
    210   //free(sld); 
    211   //free(thick_inter); 
    212   //free(thick); 
    213   //free(fun_coef); 
    214164 
    215165  return (f2); 
    … …  
    222172 * @return: function value 
    223173 */ 
    224 double Iq(double q, 
     174static double Iq(double q, 
    225175    int n_shells, 
    226     double thick_inter[], 
    227     double func_inter[], 
     176    int npts_inter, 
     177    double radius_core, 
    228178    double sld_core, 
    229179    double sld_solvent, 
    230180    double sld_flat[], 
    231181    double thick_flat[], 
    232     double nu_inter[], 
    233     int npts_inter, 
    234     double core_radius 
    235     ) { 
     182    double func_inter[], 
     183    double thick_inter[], 
     184    double nu_inter[] ) { 
    236185 
    237186    //printf("Number of points %d\n",npts_inter); 
    238187    double intensity; 
    239     //TODO: Remove this container at later stage. It is only kept to minimize stupid errors now 
     188    //TODO: Remove this container at later stage. 
    240189    double dp[60]; 
    241190    dp[0] = n_shells; 
    242191    //This is scale will also have to be removed at some stage 
    243192    dp[1] = 1.0; 
    244     dp[2] = thick_inter_0; 
    245     dp[3] = func_inter_0; 
     193    dp[2] = thick_inter[0]; 
     194    dp[3] = func_inter[0]; 
    246195    dp[4] = sld_core; 
    247196    dp[5] = sld_solvent; 
    248197    dp[6] = 0.0; 
    249198 
    250     for (i=0; i<n; i++){ 
     199    for (int i=1; i<n_shells; i++){ 
    251200        dp[i+7] = sld_flat[i]; 
    252201        dp[i+17] = thick_inter[i]; 
    … …  
    257206 
    258207    dp[57] = npts_inter; 
    259     dp[58] = nu_inter_0; 
    260     dp[59] = rad_core_0; 
     208    dp[58] = nu_inter[0]; 
     209    dp[59] = radius_core; 
    261210 
    262211    intensity = 1.0e-4*sphere_sld_kernel(dp,q); 
    … …  
    271220 * @return: function value 
    272221 */ 
    273 double Iqxy(double qx, double qy, 
     222 
     223/*static double Iqxy(double qx, double qy, 
    274224    int n_shells, 
    275     double thick_inter[], 
    276     double func_inter[], 
     225    int npts_inter, 
     226    double radius_core 
    277227    double sld_core, 
    278228    double sld_solvent, 
    279229    double sld_flat[], 
    280230    double thick_flat[], 
     231    double func_inter[], 
     232    double thick_inter[], 
    281233    double nu_inter[], 
    282     int npts_inter, 
    283     double core_radius 
    284234    ) { 
    285235 
    286236    double q = sqrt(qx*qx + qy*qy); 
    287     return Iq(q, n_shells, thick_inter_0, func_inter_0, core0_sld, solvent_sld, 
    288     flat1_sld, flat2_sld, flat3_sld, flat4_sld, flat5_sld, 
    289     flat6_sld, flat7_sld, flat8_sld, flat9_sld, flat10_sld, 
    290     thick_inter_1, thick_inter_2, thick_inter_3, thick_inter_4, thick_inter_5, 
    291     thick_inter_6, thick_inter_7, thick_inter_8, thick_inter_9, thick_inter_10, 
    292     thick_flat_1, thick_flat_2, thick_flat_3, thick_flat_4, thick_flat_5, 
    293     thick_flat_6, thick_flat_7, thick_flat_8, thick_flat_9, thick_flat_10, 
    294     func_inter_1, func_inter_2, func_inter_3, func_inter_4, func_inter_5, 
    295     func_inter_6, func_inter_7, func_inter_8, func_inter_9, func_inter_10, 
    296     nu_inter_1, nu_inter_2, nu_inter_3, nu_inter_4, nu_inter_5, 
    297     nu_inter_6, nu_inter_7, nu_inter_8, nu_inter_9, nu_inter_10, 
    298     npts_inter, nu_inter_0, rad_core_0); 
    299  
    300     //TODO: Check if evalute rphi is not needed? 
    301  
    302 } 
    303  
     237    return Iq(q, n_shells, npts_inter, radius_core, sld_core, sld_solvent, 
     238    sld_flat[], thick_flat[], func_inter[], thick_inter[], nu_inter[]) 
     239}*/ 
     240 

  • sasmodels/models/spherical_sld.py

    rd2bb604 r4605bf10  
    170170# pylint: disable=bad-whitespace, line-too-long 
    171171#            ["name", "units", default, [lower, upper], "type", "description"], 
    172 parameters = [["n",                "",               1,      [0, 9],         "", "number of shells"], 
     172parameters = [["n_shells",                "",               1,      [0, 9],         "", "number of shells"], 
     173              ["npts_inter",       "",               35,     [0, 35],        "", "number of points in each sublayer Must be odd number"], 
    173174              ["radius_core",      "Ang",            50.0,   [0, inf],       "", "intern layer thickness"], 
    174175              ["sld_core",         "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"], 
    175               ["sld_flat[n]",      "1e-6/Ang^2",     4.06,   [-inf, inf],    "", "sld function flat"], 
    176               ["thick_flat[n]",    "Ang",            100.0,  [0, inf],       "", "flat layer_thickness"], 
    177               ["func_inter[n]",    "",               0,      [0, 4],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"], 
    178               ["thick_inter[n]",   "Ang",            50.0,   [0, inf],       "", "intern layer thickness"], 
    179               ["inter_nu[n]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"], 
    180               ["npts_inter",       "",               35,     [0, 35],        "", "number of points in each sublayer Must be odd number"], 
    181176              ["sld_solvent",      "1e-6/Ang^2",     1.0,    [-inf, inf],    "", "sld function solvent"], 
     177              ["sld_flat[n_shells]",      "1e-6/Ang^2",     4.06,   [-inf, inf],    "", "sld function flat"], 
     178              ["thick_flat[n_shells]",    "Ang",            100.0,  [0, inf],       "", "flat layer_thickness"], 
     179              ["func_inter[n_shells]",    "",               0,      [0, 4],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"], 
     180              ["thick_inter[n_shells]",   "Ang",            50.0,   [0, inf],       "", "intern layer thickness"], 
     181              ["nu_inter[n_shells]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"], 
    182182              ] 
    183183# pylint: enable=bad-whitespace, line-too-long 
    184 #source = ["lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"] 
    185  
    186 def Iq(q, *args, **kw): 
    187     return q 
    188  
    189 def Iqxy(qx, *args, **kw): 
    190     return qx 
     184source = ["lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"] 
     185 
     186#def Iq(q, *args, **kw): 
     187#    return q 
     188 
     189#def Iqxy(qx, *args, **kw): 
     190#    return qx 
     191 
     192def profile(n_shells, radius_core,  sld_core,  sld_solvent, sld_flat, 
     193            thick_flat, func_inter, thick_inter, nu_inter, npts_inter): 
     194    """ 
     195    Returns shape profile with x=radius, y=SLD. 
     196    """ 
     197 
     198    z = [] 
     199    beta = [] 
     200    z0 = 0 
     201    # two sld points for core 
     202    z.append(0) 
     203    beta.append(sld_core) 
     204    z.append(radius_core) 
     205    beta.append(sld_core) 
     206    z0 += radius_core 
     207 
     208    for i in range(1, n_shells+2): 
     209        dz = thick_inter[i-1]/npts_inter 
     210        # j=0 for interface, j=1 for flat layer 
     211        for j in range(0, 2): 
     212            # interation for sub-layers 
     213            for n_s in range(0, npts_inter+1): 
     214                if j == 1: 
     215                    if i == n_shells+1: 
     216                        break 
     217                    # shift half sub thickness for the first point 
     218                    z0 -= dz#/2.0 
     219                    z.append(z0) 
     220                    #z0 -= dz/2.0 
     221                    z0 += thick_flat[i] 
     222                    sld_i = sld_flat[i] 
     223                    beta.append(sld_flat[i]) 
     224                    dz = 0 
     225                else: 
     226                    nu = nu_inter[i-1] 
     227                    # decide which sld is which, sld_r or sld_l 
     228                    if i == 1: 
     229                        sld_l = sld_core 
     230                    else: 
     231                        sld_l = sld_flat[i-1] 
     232                    if i == n_shells+1: 
     233                        sld_r = sld_solvent 
     234                    else: 
     235                        sld_r = sld_flat[i] 
     236                    # get function type 
     237                    func_idx = func_inter[i-1] 
     238                    # calculate the sld 
     239                    sld_i = intersldfunc(func_idx, npts_inter, n_s, nu, 
     240                                            sld_l, sld_r) 
     241                # append to the list 
     242                z.append(z0) 
     243                beta.append(sld_i) 
     244                z0 += dz 
     245                if j == 1: 
     246                    break 
     247    z.append(z0) 
     248    beta.append(sld_solvent) 
     249    z_ext = z0/5.0 
     250    z.append(z0+z_ext) 
     251    beta.append(sld_solvent) 
     252    # return sld profile (r, beta) 
     253    return np.asarray(z), np.asarray(beta)*1e-6 
     254 
     255def ER(core_radius, n, thickness): 
     256    return np.sum(thickness[:n[0]], axis=0) + core_radius 
     257 
     258def VR(core_radius, n, thickness): 
     259    return 1.0, 1.0 
    191260 
    192261 
    … …  
    194263    n=4, 
    195264    scale=1.0, 
    196     solvent_sld=1.0, 
     265    sld_solvent=1.0, 
    197266    background=0.0, 
    198267    npts_inter=35.0, 
Note: See TracChangeset for help on using the changeset viewer.