Changes in / [02e70ff:6ddd6e0] in sasmodels


Ignore:
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • doc/genmodel.py

    raa2edb2 rc094758  
    1010# Convert ../sasmodels/models/name.py to name 
    1111model_name = os.path.basename(sys.argv[1])[:-3] 
     12model_info = core.load_model_info(model_name) 
     13model = core.build_model(model_info) 
    1214 
    1315# Load the doc string from the module definition file and store it in rst 
    14 docstr = generate.make_doc(core.load_model_info(model_name)) 
    15      
    16 # Generate automatically plot of the model and add it to rst documentation 
     16docstr = generate.make_doc(model_info) 
    1717 
    18 info = core.load_model_info(model_name) 
    1918 
    2019# Calculate 1D curve for default parameters 
    21 pars = dict((p[0], p[2]) for p in info['parameters']) 
     20pars = dict((p.name, p.default) for p in model_info['parameters']) 
    2221 
    2322# Plotting ranges and options 
    2423opts = { 
    25         'xscale'    : 'log', 
    26         'yscale'    : 'log' if not info['structure_factor'] else 'linear', 
    27         'qmin'      : 0.001, 
    28         'qmax'      : 1.0, 
    29         'nq'        : 1000, 
    30         'nq2d'      : 100, 
     24    'xscale'    : 'log', 
     25    'yscale'    : 'log' if not model_info['structure_factor'] else 'linear', 
     26    'zscale'    : 'log' if not model_info['structure_factor'] else 'linear', 
     27    'q_min'     : 0.001, 
     28    'q_max'     : 1.0, 
     29    'nq'        : 1000, 
     30    'nq2d'      : 400, 
     31    'vmin'      : 1e-3,  # floor for the 2D data results 
     32    'qx_max'    : 0.5, 
    3133} 
    3234 
    33 qmin, qmax, nq = opts['qmin'], opts['qmax'], opts['nq'] 
    34 qmin = math.log10(qmin) 
    35 qmax = math.log10(qmax) 
    36 q = np.logspace(qmin, qmax, nq) 
    37 data = empty_data1D(q) 
    38 model = core.load_model(model_name) 
    39 calculator = DirectModel(data, model) 
    40 Iq1D = calculator() 
    4135 
    42 # TO DO: Generation of 2D plots  
    43 # Problem in sasmodels.direct_model._calc_theory 
    44 # There self._kernel.q_input.nq  gets a value of 0 in the 2D case 
    45 # and returns a 0 numpy array (it does not call the C code) 
     36def plot_1d(model, opts, ax): 
     37    q_min, q_max, nq = opts['q_min'], opts['q_max'], opts['nq'] 
     38    q_min = math.log10(q_min) 
     39    q_max = math.log10(q_max) 
     40    q = np.logspace(q_min, q_max, nq) 
     41    data = empty_data1D(q) 
     42    calculator = DirectModel(data, model) 
     43    Iq1D = calculator() 
    4644 
    47 # If 2D model, compute 2D image 
    48 #if info['has_2d'] != []: 
    49 #    qmax, nq2d = opts['qmax'], opts['nq2d'] 
    50 #    data2d = empty_data2D(np.linspace(-qmax, qmax, nq2d), resolution=0.0)  
    51 #    #model = core.load_model(model_name) 
    52 #    calculator = DirectModel(data2d, model) 
    53 #    Iq2D = calculator() 
     45    ax.plot(q, Iq1D, color='blue', lw=2, label=model_info['name']) 
     46    ax.set_xlabel(r'$Q \/(\AA^{-1})$') 
     47    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$') 
     48    ax.set_xscale(opts['xscale']) 
     49    ax.set_yscale(opts['yscale']) 
     50    #ax.legend(loc='best') 
     51 
     52def plot_2d(model, opts, ax): 
     53    qx_max, nq2d = opts['qx_max'], opts['nq2d'] 
     54    q = np.linspace(-qx_max, qx_max, nq2d) 
     55    data2d = empty_data2D(q, resolution=0.0) 
     56    calculator = DirectModel(data2d, model) 
     57    Iq2D = calculator() #background=0) 
     58    Iq2D = Iq2D.reshape(nq2d, nq2d) 
     59    if opts['zscale'] == 'log': 
     60        Iq2D = np.log(np.clip(Iq2D, opts['vmin'], np.inf)) 
     61    h = ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='upper', 
     62           extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=ice_cm()) 
     63    # , vmin=vmin, vmax=vmax) 
     64    ax.set_xlabel(r'$Q_x \/(\AA^{-1})$') 
     65    ax.set_ylabel(r'$Q_y \/(\AA^{-1})$') 
     66 
     67def ice_cm(): 
     68    from matplotlib._cm import _Blues_data 
     69    from matplotlib import colors 
     70    from matplotlib import rcParams 
     71    def from_white(segments): 
     72        scale = 1.0/segments[0][1] 
     73        return [(k, v*scale, w*scale) for k, v, w in segments] 
     74    ice_data = dict((k,from_white(v)) for k,v in _Blues_data.items()) 
     75    ice = colors.LinearSegmentedColormap("ice", ice_data, rcParams['image.lut']) 
     76    return ice 
     77 
    5478 
    5579# Generate image (comment IF for 1D/2D for the moment) and generate only 1D 
    56 #if info['has_2d'] == []: 
    57 #    fig = plt.figure() 
    58 #    ax = fig.add_subplot(1,1,1) 
    59 #    ax.plot(q, Iq1D, color='blue', lw=2, label=model_name) 
    60 #    ax.set_xlabel(r'$Q \/(\AA^{-1})$') 
    61 #    ax.set_xscale(opts['xscale'])    
    62 #    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$') 
    63 #    ax.set_yscale(opts['yscale'])   
    64 #    ax.legend() 
    65 #else: 
    66 #    # need figure with 1D + 2D 
    67 #    pass 
    68 fig = plt.figure() 
    69 ax = fig.add_subplot(1,1,1) 
    70 ax.plot(q, Iq1D, color='blue', lw=2, label=info['name']) 
    71 ax.set_xlabel(r'$Q \/(\AA^{-1})$') 
    72 ax.set_xscale(opts['xscale'])    
    73 ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$') 
    74 ax.set_yscale(opts['yscale'])   
    75 ax.legend() 
    76   
     80fig_height = 3.0 # in 
     81fig_left = 0.6 # in 
     82fig_right = 0.5 # in 
     83fig_top = 0.6*0.25 # in 
     84fig_bottom = 0.6*0.75 
     85if model_info['has_2d']: 
     86    plot_height = fig_height - (fig_top+fig_bottom) 
     87    plot_width = plot_height 
     88    fig_width = 2*(plot_width + fig_left + fig_right) 
     89    aspect = (fig_width, fig_height) 
     90    ratio = aspect[0]/aspect[1] 
     91    ax_left = fig_left/fig_width 
     92    ax_bottom = fig_bottom/fig_height 
     93    ax_height = plot_height/fig_height 
     94    ax_width = ax_height/ratio # square axes 
     95    fig = plt.figure(figsize=aspect) 
     96    ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height]) 
     97    plot_2d(model, opts, ax2d) 
     98    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
     99    #ax.set_aspect('square') 
     100else: 
     101    plot_height = fig_height - (fig_top+fig_bottom) 
     102    plot_width = (1+np.sqrt(5))/2*fig_height 
     103    fig_width = plot_width + fig_left + fig_right 
     104    ax_left = fig_left/fig_width 
     105    ax_bottom = fig_bottom/fig_height 
     106    ax_width = plot_width/fig_width 
     107    ax_height = plot_height/fig_height 
     108    aspect = (fig_width, fig_height) 
     109    fig = plt.figure(figsize=aspect) 
     110    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
     111plot_1d(model, opts, ax1d) 
    77112 
    78113# Save image in model/img 
    79114figname = model_name + '_autogenfig.png' 
    80115filename = os.path.join('model', 'img', figname) 
    81 plt.savefig(filename) 
     116plt.savefig(filename, bbox_inches='tight') 
     117#print "figure saved in",filename 
    82118 
    83119# Auto caption for figure 
    84120captionstr = '\n' 
    85 captionstr += '.. figure:: img/' + model_name + '_autogenfig.png\n' 
     121captionstr += '.. figure:: img/' + model_info['id'] + '_autogenfig.png\n' 
    86122captionstr += '\n' 
    87 #if info['has_2d'] == []: 
    88 #    captionstr += '    1D plot corresponding to the default parameters of the model.\n' 
    89 #else: 
    90 #    captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n' 
    91123captionstr += '    1D plot corresponding to the default parameters of the model.\n' 
    92124captionstr += '\n' 
    … …  
    102134else: 
    103135    print '------------------------------------------------------------------' 
    104     print 'References NOT FOUND for model: ', model_name 
     136    print 'References NOT FOUND for model: ', model_info['id'] 
    105137    print '------------------------------------------------------------------' 
    106138    docstr = docstr + captionstr 

  • explore/J1c.py

    r0a6da3c rcbd37a7  
    99 
    1010 
    11 SHOW_DIFF = True  # True if show diff rather than function value 
     11SHOW_DIFF = True # True if show diff rather than function value 
     12#SHOW_DIFF = False # True if show diff rather than function value 
    1213LINEAR_X = False  # True if q is linearly spaced instead of log spaced 
     14#LINEAR_X = True # True if q is linearly spaced instead of log spaced 
     15FUNCTION = "2*J1(x)/x" 
    1316 
    14 def mp_J1c(vec, bits=500): 
     17def mp_fn(vec, bits=500): 
    1518    """ 
    1619    Direct calculation using sympy multiprecision library. 
    1720    """ 
    1821    with mp.workprec(bits): 
    19         return [_mp_J1c(mp.mpf(x)) for x in vec] 
     22        return [_mp_fn(mp.mpf(x)) for x in vec] 
    2023 
    21 def _mp_J1c(x): 
     24def _mp_fn(x): 
    2225    """ 
    23     Helper funciton for mp_j1c 
     26    Actual function that gets evaluated.  The caller just vectorizes. 
    2427    """ 
    2528    return mp.mpf(2)*mp.j1(x)/x 
    2629 
    27 def np_J1c(x, dtype): 
     30def np_fn(x, dtype): 
    2831    """ 
    2932    Direct calculation using scipy. 
    … …  
    3336    return np.asarray(2, dtype)*J1(x)/x 
    3437 
    35 def cephes_J1c(x, dtype, n): 
     38def sasmodels_fn(x, dtype, platform='ocl'): 
    3639    """ 
    3740    Calculation using pade approximant. 
    3841    """ 
    39     f = np.float64 if np.dtype(dtype) == np.float64 else np.float32 
    40     x = np.asarray(x, dtype) 
    41     ans = np.empty_like(x) 
    42     ax = abs(x) 
    43  
    44     # Branch a 
    45     a_idx = ax < f(8.0) 
    46     a_xsq = x[a_idx]**2 
    47     a_coeff1 = list(reversed((72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606))) 
    48     a_coeff2 = list(reversed((144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0))) 
    49     a_ans1 = np.polyval(np.asarray(a_coeff1[n:], dtype), a_xsq) 
    50     a_ans2 = np.polyval(np.asarray(a_coeff2[n:], dtype), a_xsq) 
    51     ans[a_idx] = f(2.0)*a_ans1/a_ans2 
    52  
    53     # Branch b 
    54     b_idx = ~a_idx 
    55     b_ax = ax[b_idx] 
    56     b_x = x[b_idx] 
    57  
    58     b_y = f(64.0)/(b_ax**2) 
    59     b_xx = b_ax - f(2.356194491) 
    60     b_coeff1 = list(reversed((1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6))) 
    61     b_coeff2 = list(reversed((0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6))) 
    62     b_ans1 = np.polyval(np.asarray(b_coeff1[n:], dtype),b_y) 
    63     b_ans2 = np.polyval(np.asarray(b_coeff2[n:], dtype), b_y) 
    64     b_sn, b_cn = np.sin(b_xx), np.cos(b_xx) 
    65     ans[b_idx] = np.sign(b_x)*np.sqrt(f(0.636619772)/b_ax) * (b_cn*b_ans1 - (f(8.0)/b_ax)*b_sn*b_ans2)*f(2.0)/b_x 
    66  
    67     return ans 
    68  
    69 def div_J1c(x, dtype): 
    70     f = np.float64 if np.dtype(dtype) == np.float64 else np.float32 
    71     x = np.asarray(x, dtype) 
    72     return f(2.0)*np.asarray([_J1(xi, f)/xi for xi in x], dtype) 
    73  
    74 def _J1(x, f): 
    75     ax = abs(x) 
    76     if ax < f(8.0): 
    77         y = x*x 
    78         ans1 = x*(f(72362614232.0) 
    79                   + y*(f(-7895059235.0) 
    80                   + y*(f(242396853.1) 
    81                   + y*(f(-2972611.439) 
    82                   + y*(f(15704.48260) 
    83                   + y*(f(-30.16036606))))))) 
    84         ans2 = (f(144725228442.0) 
    85                   + y*(f(2300535178.0) 
    86                   + y*(f(18583304.74) 
    87                   + y*(f(99447.43394) 
    88                   + y*(f(376.9991397) 
    89                   + y))))) 
    90         return ans1/ans2 
    91     else: 
    92         y = f(64.0)/(ax*ax) 
    93         xx = ax - f(2.356194491) 
    94         ans1 = (f(1.0) 
    95                   + y*(f(0.183105e-2) 
    96                   + y*(f(-0.3516396496e-4) 
    97                   + y*(f(0.2457520174e-5) 
    98                   + y*f(-0.240337019e-6))))) 
    99         ans2 = (f(0.04687499995) 
    100                   + y*(f(-0.2002690873e-3) 
    101                   + y*(f(0.8449199096e-5) 
    102                   + y*(f(-0.88228987e-6) 
    103                   + y*f(0.105787412e-6))))) 
    104         sn, cn = np.sin(xx), np.cos(xx) 
    105         ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) 
    106         return -ans if (x < f(0.0)) else ans 
     42    from sasmodels import core, data, direct_model 
     43    model = core.load_model('bessel', dtype=dtype) 
     44    calculator = direct_model.DirectModel(data.empty_data1D(x), model) 
     45    return calculator(background=0) 
    10746 
    10847def plotdiff(x, target, actual, label): 
    … …  
    11352    """ 
    11453    if SHOW_DIFF: 
    115         err = np.clip(abs((target-actual)/target), 0, 1) 
     54        err = abs((target-actual)/target) 
     55        #err = np.clip(err, 0, 1) 
    11656        pylab.loglog(x, err, '-', label=label) 
    11757    else: 
    … …  
    11959        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label) 
    12060 
    121 def compare(x, precision): 
     61def compare(x, precision, target): 
    12262    r""" 
    12363    Compare the different computation methods using the given precision. 
    12464    """ 
    125     target = np.asarray(mp_J1c(x), 'double') 
    126     #plotdiff(x, target, mp_J1c(x, 11), 'mp 11 bits') 
    127     plotdiff(x, target, np_J1c(x, precision), 'direct '+precision) 
    128     plotdiff(x, target, cephes_J1c(x, precision, 0), 'cephes '+precision) 
    129     #plotdiff(x, target, cephes_J1c(x, precision, 1), 'cephes '+precision) 
    130     #plotdiff(x, target, div_J1c(x, precision), 'cephes 2 J1(x)/x '+precision) 
     65    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits') 
     66    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision) 
     67    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision) 
    13168    pylab.xlabel("qr (1/Ang)") 
    13269    if SHOW_DIFF: 
    13370        pylab.ylabel("relative error") 
    13471    else: 
    135         pylab.ylabel("2 J1(x)/x") 
     72        pylab.ylabel(FUNCTION) 
    13673        pylab.semilogx(x, target,  '-', label="true value") 
    13774    if LINEAR_X: 
    … …  
    14784    else: 
    14885        qr = np.logspace(-3,5,400) 
     86    target = np.asarray(mp_fn(qr), 'double') 
    14987    pylab.subplot(121) 
    150     compare(qr, 'single') 
     88    compare(qr, 'single', target) 
    15189    pylab.legend(loc='best') 
    15290    pylab.subplot(122) 
    153     compare(qr, 'double') 
     91    compare(qr, 'double', target) 
    15492    pylab.legend(loc='best') 
    155     pylab.suptitle('2 J1(x)/x') 
     93    pylab.suptitle(FUNCTION) 
    15694 
    15795if __name__ == "__main__": 

  • sasmodels/data.py

    r84db7a5 rc094758  
    178178        self.qy_data, self.dqy_data = y, dy 
    179179        self.data, self.err_data = z, dz 
    180         self.mask = (~np.isnan(z) if z is not None 
    181                      else np.ones_like(x) if x is not None 
     180        self.mask = (np.isnan(z) if z is not None 
     181                     else np.zeros_like(x, dtype='bool') if x is not None 
    182182                     else None) 
    183183        self.q_data = np.sqrt(x**2 + y**2) 

  • sasmodels/kernelcl.py

    r17bbadd rc094758  
    367367        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 
    368368        context = env.get_context(self.dtype) 
     369        self.global_size = [self.q_vectors[0].size] 
     370        #print("creating inputs of size", self.global_size) 
    369371        self.q_buffers = [ 
    370372            cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 
    371373            for q in self.q_vectors 
    372374        ] 
    373         self.global_size = [self.q_vectors[0].size] 
    374375 
    375376    def release(self): 

  • sasmodels/models/bessel.py

    r07142f3 rcbd37a7  
    6767#Bessel 
    6868parameters = [ 
     69    ["ignored", "", 0.0, [-inf, inf], "", "no parameterless functions"], 
    6970             ] 
    7071 
    71 source = ["lib/polevl.c", "lib/j1d.c"] 
     72source = ["lib/polevl.c", "lib/j1_cephes.c"] 
    7273 
    7374# No volume normalization despite having a volume parameter 
    … …  
    7778 
    7879Iq = """ 
    79     return j1(q); 
     80    return 2.0*j1(q)/q; 
    8081    """ 
    8182 

  • sasmodels/models/lib/j0_cephes.c

    rbfef528 r094e320  
    4444 */ 
    4545 
    46 /*                                                      y0.c 
    47  * 
    48  *      Bessel function of the second kind, order zero 
    49  * 
    50  * 
    51  * 
    52  * SYNOPSIS: 
    53  * 
    54  * double x, y, y0(); 
    55  * 
    56  * y = y0( x ); 
    57  * 
    58  * 
    59  * 
    60  * DESCRIPTION: 
    61  * 
    62  * Returns Bessel function of the second kind, of order 
    63  * zero, of the argument. 
    64  * 
    65  * The domain is divided into the intervals [0, 5] and 
    66  * (5, infinity). In the first interval a rational approximation 
    67  * R(x) is employed to compute 
    68  *   y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI. 
    69  * Thus a call to j0() is required. 
    70  * 
    71  * In the second interval, the Hankel asymptotic expansion 
    72  * is employed with two rational functions of degree 6/6 
    73  * and 7/7. 
    74  * 
    75  * 
    76  * 
    77  * ACCURACY: 
    78  * 
    79  *  Absolute error, when y0(x) < 1; else relative error: 
    80  * 
    81  * arithmetic   domain     # trials      peak         rms 
    82  *    DEC       0, 30        9400       7.0e-17     7.9e-18 
    83  *    IEEE      0, 30       30000       1.3e-15     1.6e-16 
    84  * 
    85  */ 
    86  
    8746 
    8847/* 
    … …  
    9554 
    9655double j0( double ); 
    97  
    9856double j0(double x) { 
    9957 
    … …  
    291249 
    292250    q = 1.0/x; 
    293     w = sqrtf(q); 
     251    w = sqrt(q); 
    294252 
    295253    p = w * polevl( q, MO, 7); 
    296254    w = q*q; 
    297255    xn = q * polevl( w, PH, 7) - PIO4F; 
    298     p = p * cosf(xn + xx); 
     256    p = p * cos(xn + xx); 
    299257    return(p); 
    300258#endif 

  • sasmodels/models/lib/j1_cephes.c

    rbfef528 re2af2a9  
    3232 *    IEEE      0, 30       30000       2.6e-16     1.1e-16 
    3333 * 
    34  * 
    35  */ 
    36 /*                                                      y1.c 
    37  * 
    38  *      Bessel function of second kind of order one 
    39  * 
    40  * 
    41  * 
    42  * SYNOPSIS: 
    43  * 
    44  * double x, y, y1(); 
    45  * 
    46  * y = y1( x ); 
    47  * 
    48  * 
    49  * 
    50  * DESCRIPTION: 
    51  * 
    52  * Returns Bessel function of the second kind of order one 
    53  * of the argument. 
    54  * 
    55  * The domain is divided into the intervals [0, 8] and 
    56  * (8, infinity). In the first interval a 25 term Chebyshev 
    57  * expansion is used, and a call to j1() is required. 
    58  * In the second, the asymptotic trigonometric representation 
    59  * is employed using two rational functions of degree 5/5. 
    60  * 
    61  * 
    62  * 
    63  * ACCURACY: 
    64  * 
    65  *                      Absolute error: 
    66  * arithmetic   domain      # trials      peak         rms 
    67  *    DEC       0, 30       10000       8.6e-17     1.3e-17 
    68  *    IEEE      0, 30       30000       1.0e-15     1.3e-16 
    69  * 
    70  * (error criterion relative when |y1| > 1). 
    7134 * 
    7235 */ 

  • sasmodels/models/lib/polevl.c

    r3936ad3 re2af2a9  
    5050Direct inquiries to 30 Frost Street, Cambridge, MA 02140 
    5151*/ 
     52 
    5253double polevl( double x, double coef[8], int N ); 
    5354double p1evl( double x, double coef[8], int N ); 
Note: See TracChangeset for help on using the changeset viewer.