Changeset fdc538a in sasmodels


Ignore:
Timestamp:
Jan 27, 2016 5:47:17 PM (9 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:
09ebe8c
Parents:
d15a908
Message:

delint

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/resolution.py

    r5258859 rfdc538a  
    135135    """ 
    136136    #print("apply shapes", theory.shape, weight_matrix.shape) 
    137     Iq = np.dot(theory[None,:], weight_matrix) 
     137    Iq = np.dot(theory[None, :], weight_matrix) 
    138138    #print("result shape",Iq.shape) 
    139139    return Iq.flatten() 
     
    153153    # neither trapezoid nor Simpson's rule improved the accuracy. 
    154154    edges = bin_edges(q_calc) 
    155     edges[edges<0.0] = 0.0 # clip edges below zero 
    156     G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:] ) 
     155    edges[edges < 0.0] = 0.0 # clip edges below zero 
     156    G = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    157157    weights = G[1:] - G[:-1] 
    158     weights /= np.sum(weights, axis=0)[None,:] 
     158    weights /= np.sum(weights, axis=0)[None, :] 
    159159    return weights 
    160160 
     
    287287    # The current algorithm is a midpoint rectangle rule. 
    288288    q_edges = bin_edges(q_calc) # Note: requires q > 0 
    289     q_edges[q_edges<0.0] = 0.0 # clip edges below zero 
     289    q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 
    290290    weights = np.zeros((len(q), len(q_calc)), 'd') 
    291291 
     
    306306            #print(qi - h, qi + h) 
    307307            #print(in_x + abs_x) 
    308             weights[i,:] = (in_x + abs_x) * np.diff(q_edges) / (2*h) 
     308            weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h) 
    309309        else: 
    310310            L = n_height 
    311311            for k in range(-L, L+1): 
    312                 weights[i,:] += _q_perp_weights(q_edges, qi+k*h/L, w) 
    313             weights[i,:] /= 2*L + 1 
     312                weights[i, :] += _q_perp_weights(q_edges, qi+k*h/L, w) 
     313            weights[i, :] /= 2*L + 1 
    314314 
    315315    return weights.T 
     
    358358    log-scaled data. 
    359359    """ 
    360     if len(x) < 2 or (np.diff(x)<0).any(): 
     360    if len(x) < 2 or (np.diff(x) < 0).any(): 
    361361        raise ValueError("Expected bins to be an increasing set") 
    362362    edges = np.hstack([ 
     
    373373    """ 
    374374    step = np.diff(q) 
    375     index = step>max_step 
     375    index = step > max_step 
    376376    if np.any(index): 
    377377        inserts = [] 
    378         for q_i,step_i in zip(q[:-1][index],step[index]): 
     378        for q_i, step_i in zip(q[:-1][index], step[index]): 
    379379            n = np.ceil(step_i/max_step) 
    380             inserts.extend(q_i + np.arange(1,n)*(step_i/n)) 
     380            inserts.extend(q_i + np.arange(1, n)*(step_i/n)) 
    381381        # Extend a couple of fringes beyond the end of the data 
    382         inserts.extend(q[-1] + np.arange(1,8)*max_step) 
    383         q_calc = np.sort(np.hstack((q,inserts))) 
     382        inserts.extend(q[-1] + np.arange(1, 8)*max_step) 
     383        q_calc = np.sort(np.hstack((q, inserts))) 
    384384    else: 
    385385        q_calc = q 
     
    399399    if q_min < q[0]: 
    400400        if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    401         n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1]>q[0] else 15 
     401        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15 
    402402        q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 
    403403    else: 
    404404        q_low = [] 
    405405    if q_max > q[-1]: 
    406         n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1]>q[-2] else 15 
     406        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15 
    407407        q_high = np.linspace(q[-1], q_max, n_high+1)[1:] 
    408408    else: 
     
    452452        if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    453453        n_low = log_delta_q * (log(q[0])-log(q_min)) 
    454         q_low  = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 
     454        q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 
    455455    else: 
    456456        q_low = [] 
     
    489489    that make it slow to evaluate but give it good accuracy. 
    490490    """ 
    491     from scipy.integrate import romberg, dblquad 
     491    from scipy.integrate import romberg 
    492492 
    493493    if any(k not in form.info['defaults'] for k in pars.keys()): 
     
    520520            result[i] = r/(2*h) 
    521521        else: 
    522             w_grid = np.linspace(0, w, 21)[None,:] 
    523             h_grid = np.linspace(-h, h, 23)[:,None] 
     522            w_grid = np.linspace(0, w, 21)[None, :] 
     523            h_grid = np.linspace(-h, h, 23)[:, None] 
    524524            u = sqrt((qi+h_grid)**2 + w_grid**2) 
    525525            Iu = np.interp(u, q_calc, Iq) 
    526526            #print(np.trapz(Iu, w_grid, axis=1)) 
    527             Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:,0]) 
     527            Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:, 0]) 
    528528            result[i] = Is / (2*h*w) 
    529             """ 
    530             r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, 
    531                              args=(qi,)) 
    532             result[i] = r/(w*2*h) 
    533             """ 
     529            # from scipy.integrate import dblquad 
     530            # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, 
     531            #                  args=(qi,)) 
     532            # result[i] = r/(w*2*h) 
    534533 
    535534    # r should be [float, ...], but it is [array([float]), array([float]),...] 
     
    553552 
    554553    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 
    555     r = [romberg(_fn, max(qi-nsigma*dqi,1e-10*q[0]), qi+nsigma*dqi, args=(qi, dqi), 
    556                  divmax=100, vec_func=True, tol=0, rtol=1e-8) 
    557          for qi,dqi in zip(q,q_width)] 
     554    r = [romberg(_fn, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi, 
     555                 args=(qi, dqi), divmax=100, vec_func=True, tol=0, rtol=1e-8) 
     556         for qi, dqi in zip(q, q_width)] 
    558557    return np.asarray(r).flatten() 
    559558 
     
    595594        theory = self.Iq(resolution.q_calc) 
    596595        output = resolution.apply(theory) 
    597         answer = [ 9.0618, 8.6402, 8.1187, 7.1392, 6.1528, 
    598                    5.5555, 4.5584, 3.5606, 2.5623, 2.0000 ] 
     596        answer = [ 
     597            9.0618, 8.6402, 8.1187, 7.1392, 6.1528, 
     598            5.5555, 4.5584, 3.5606, 2.5623, 2.0000, 
     599            ] 
    599600        np.testing.assert_allclose(output, answer, atol=1e-4) 
    600601 
     
    608609        theory = 1000*self.Iq(resolution.q_calc**4) 
    609610        output = resolution.apply(theory) 
    610         answer = [ 8.85785, 8.43012, 7.92687, 6.94566, 6.03660, 
    611                    5.40363, 4.40655, 3.40880, 2.41058, 2.00000 ] 
     611        answer = [ 
     612            8.85785, 8.43012, 7.92687, 6.94566, 6.03660, 
     613            5.40363, 4.40655, 3.40880, 2.41058, 2.00000, 
     614            ] 
    612615        np.testing.assert_allclose(output, answer, atol=1e-4) 
    613616 
     
    620623        theory = self.Iq(resolution.q_calc) 
    621624        output = resolution.apply(theory) 
    622         answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ] 
     625        answer = [ 
     626            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 
     627            ] 
    623628        np.testing.assert_allclose(output, answer, atol=1e-4) 
    624629 
     
    632637        theory = self.Iq(resolution.q_calc) 
    633638        output = resolution.apply(theory) 
    634         answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ] 
     639        answer = [ 
     640            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 
     641            ] 
    635642        np.testing.assert_allclose(output, answer, atol=1e-4) 
    636643 
     
    652659        theory = 12.0-1000.0*resolution.q_calc 
    653660        output = resolution.apply(theory) 
    654         answer = [ 10.44785079, 9.84991299, 8.98101708, 
    655                   7.99906585, 6.99998311, 6.00001689, 
    656                   5.00093415, 4.01898292, 3.15008701, 2.55214921] 
     661        answer = [ 
     662            10.44785079, 9.84991299, 8.98101708, 
     663            7.99906585, 6.99998311, 6.00001689, 
     664            5.00093415, 4.01898292, 3.15008701, 2.55214921, 
     665            ] 
    657666        np.testing.assert_allclose(output, answer, atol=1e-8) 
    658667 
     
    783792            } 
    784793        form = load_model('ellipsoid', dtype='double') 
    785         q = np.logspace(log10(4e-5),log10(2.5e-2), 68) 
     794        q = np.logspace(log10(4e-5), log10(2.5e-2), 68) 
    786795        width, height = 0.117, 0. 
    787796        resolution = Slit1D(q, width=width, height=height) 
     
    10711080    #h = 0.0277790 
    10721081    resolution = Slit1D(q, w, h) 
    1073     _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w,h)) 
     1082    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h)) 
    10741083 
    10751084def demo(): 
Note: See TracChangeset for help on using the changeset viewer.