Changeset 09e15be in sasmodels


Ignore:
Timestamp:
Jul 18, 2014 2:25:00 PM (10 years ago)
Author:
HMP1 <helen.park@…>
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:
79fcc40
Parents:
8a21ba3
Message:

Attempt at faster kernel for TEST,
updated fit.py,
errors in the kernels fixed

Files:
2 added
2 deleted
13 edited

Legend:

Unmodified
Added
Removed

  • Kernel-CapCyl.cpp

    r8a21ba3 r09e15be  
    5151            _ptvalue[i] = 0.0; 
    5252        } 
    53         if (size>1) { 
    54             _ptvalue[i] *= fabs(cos(thet*pi/180.0)); 
    55         } 
     53     //   if (size>1) { 
     54   //         _ptvalue[i] *= fabs(cos(thet*pi/180.0)); 
     55     //   } 
    5656    } 
    5757} 

  • Kernel-CoreShellCylinder.cpp

    r496b252 r09e15be  
    1111        real pi = 4.0*atan(1.0); 
    1212        real theta = axis_theta*pi/180.0; 
    13         real cyl_x = cos(theta)*cos(axis_phi*pi/180.0); 
    14         real cyl_y = sin(theta); 
    15         real alpha = acos(cyl_x*qx[i]/q + cyl_y*qy[i]/q); 
     13        real alpha = acos(cos(theta)*cos(axis_phi*pi/180.0)*qx[i]/q + sin(theta)*qy[i]/q); 
    1614 
    1715        if (alpha == 0.0){ 
    … …  
    3028 
    3129        if (besarg1 == 0.0){be1 = 0.5;} 
    32         else{ 
    33             be1 = NR_BessJ1(besarg1)/besarg1; 
    34         } 
     30        else{be1 = NR_BessJ1(besarg1)/besarg1;} 
     31 
    3532        if (besarg2 == 0.0){be2 = 0.5;} 
    36         else{ 
    37             be2 = NR_BessJ1(besarg2)/besarg2; 
    38         } 
    39         if (sinarg1 == 0.0){ 
    40             si1 = 1.0; 
    41         } 
    42         else{ 
    43             si1 = sin(sinarg1)/sinarg1; 
    44         } 
    45         if (sinarg2 == 0.0){ 
    46             si2 = 1.0; 
    47         } 
    48         else{ 
    49             si2 = sin(sinarg2)/sinarg2; 
    50         } 
    51         real tt = 2.0*vol2*(shell_sld-solvent_sld)*si2*be2 + 2.0*(pi*radius*radius*(length))*(core_sld-shell_sld)*si1*be1; 
     33        else{be2 = NR_BessJ1(besarg2)/besarg2;} 
     34 
     35        if (sinarg1 == 0.0){si1 = 1.0;} 
     36        else{si1 = sin(sinarg1)/sinarg1;} 
     37 
     38        if (sinarg2 == 0.0){si2 = 1.0;} 
     39        else{si2 = sin(sinarg2)/sinarg2;} 
     40 
     41        real tt = 2.0*vol2*(shell_sld-solvent_sld)*si2*be2+2.0*(pi*radius*radius*(length))*(core_sld-shell_sld)*si1*be1; 
    5242 
    5343        real answer = (tt*tt)*sin(alpha)/fabs(sin(alpha)); 
    54         real vol=pi*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
    55         answer = answer/vol*1.0e8*scale; 
     44        answer *= answer/(pi*(radius+thickness)*(radius+thickness)*(length+2.0*thickness))*1.0e8*scale; 
    5645 
    57         _ptvalue[i] = radius_weight*length_weight*thickness_weight*theta_weight*phi_weight*answer; 
    58         _ptvalue[i] *= pow(radius+thickness,2)*(length+2.0*thickness); 
    59  
    60         if (size>1) { 
    61         _ptvalue[i] *= fabs(cos(axis_theta*pi/180.0)); 
    62         } 
     46        _ptvalue[i] = radius_weight*length_weight*thickness_weight*theta_weight*phi_weight*answer*pow(radius+thickness,2)*(length+2.0*thickness); 
     47     //   if (size>1) { 
     48       // _ptvalue[i] *= fabs(cos(axis_theta*pi/180.0)); 
     49        //} 
    6350 
    6451    } 

  • Kernel-Cylinder.cpp

    r8a20be5 r09e15be  
    1111    if(i < count) 
    1212    { 
     13        real be=0.0; real si=0.0; 
    1314        real qq = sqrt(qx[i]*qx[i]+qy[i]*qy[i]); 
    1415 
    1516        real pi = 4.0*atan(1.0); 
    1617        real theta = cyl_theta*pi/180.0; 
    17         real phi = cyl_phi*pi/180.0; 
    18  
    19         real cyl_x = cos(theta)*cos(phi); 
    20         real cyl_y = sin(theta); 
    21         real cos_val = cyl_x*(qx[i]/qq) + cyl_y*(qy[i]/qq); 
     18        real cos_val = cos(theta)*cos(cyl_phi*pi/180.0)*(qx[i]/qq) + sin(theta)*(qy[i]/qq); 
    2219 
    2320        real alpha = acos(cos_val); 
    … …  
    2724        real besarg = qq*rr*sin(alpha); 
    2825        real siarg = qq*h/2*cos(alpha); 
    29         real be=0.0; real si=0.0; 
    3026 
    3127        real bj = NR_BessJ1(besarg); 
    32  
    3328        real d1 = qq*rr*sin(alpha); 
    3429 
    … …  
    4136        if(siarg == 0.0){ 
    4237            si = 1.0; 
    43         } 
    44         else{ 
     38        }else{ 
    4539            si = sin(siarg)*sin(siarg)/(siarg*siarg); 
    4640        } 
    47  
    4841        real form = be*si/sin(alpha); 
    4942        real answer = sub*sub*form*acos(-1.0)*rr*rr*h*1.0e8*scale; 
    5043 
    5144        _ptvalue[i] = radius_weight*length_weight*theta_weight*phi_weight*answer*pow(rr,2)*h; 
    52         if (size>1) { 
    53             _ptvalue[i] *= fabs(cos(cyl_theta*pi/180.0)); 
    54         } 
     45       // if (size>1) { 
     46         //   _ptvalue[i] *= fabs(cos(cyl_theta*pi/180.0)); 
     47      //  } 
    5548    } 
    5649} 

  • Kernel-Ellipse.cpp

    r8a21ba3 r09e15be  
    99         real pi = 4.0*atan(1.0); 
    1010         real theta = axis_theta*pi/180.0; 
    11          real cyl_x = cos(theta)*cos(axis_phi*pi/180.0); 
    12          real cyl_y = sin(theta); 
    13          real cos_val = cyl_x*(qx[i]/q) + cyl_y*(qy[i]/q); 
     11         real cos_val = cos(theta)*cos(axis_phi*pi/180.0)*(qx[i]/q) + sin(theta)*(qy[i]/q); 
    1412 
    1513         real arg = q*radius_b*sqrt(1.0+(cos_val*cos_val*(((radius_a*radius_a/(radius_b*radius_b))-1.0)))); 
    … …  
    2321 
    2422         _ptvalue[i] = radius_a_weight*radius_b_weight*axis_theta_weight*radius_a*axis_phi_weight*ret*pow(radius_b, 2); 
    25          if(size > 1){ 
    26             _ptvalue[i] *= fabs(cos(axis_theta*pi/180.0)); 
    27          } 
     23         //if(size > 1){ 
     24          //  _ptvalue[i] *= fabs(cos(axis_theta*pi/180.0)); 
     25         //} 
    2826     } 
    2927} 

  • Kernel-TriaxialEllipse.cpp

    rbe5d7df r09e15be  
    3333 
    3434        _ptvalue[i] = axisA_weight*axisB_weight*axisC_weight*theta_weight*phi_weight*psi_weight*answer*axisA*axisB*axisC; 
    35         if (size>1) 
    36         { 
    37             _ptvalue[i] *= fabs(cos(theta*pi/180.0)); 
    38         } 
     35       // if (size>1) 
     36         //   _ptvalue[i] *= fabs(cos(theta*pi/180.0)); 
     37 
    3938    } 
    4039} 

  • code_coreshellcyl.py

    r496b252 r09e15be  
    5858 
    5959        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64 
    60         for i in xrange(len(radius.weight)): 
    61             for j in xrange(len(length.weight)): 
    62                 for k in xrange(len(axis_theta.weight)): 
    63                     for l in xrange(len(axis_phi.weight)): 
    64                         for f in xrange(len(thickness.weight)): 
     60        for r in xrange(len(radius.weight)): 
     61            for l in xrange(len(length.weight)): 
     62                for at in xrange(len(axis_theta.weight)): 
     63                    for p in xrange(len(axis_phi.weight)): 
     64                        for th in xrange(len(thickness.weight)): 
    6565 
    6666                            self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, 
    67                                     real(axis_theta.value[k]), real(axis_phi.value[l]), real(thickness.value[f]), 
    68                                     real(length.value[j]), real(radius.value[i]), real(pars['scale']), 
    69                                     real(radius.weight[i]), real(length.weight[j]), real(thickness.weight[f]), 
    70                                     real(axis_theta.weight[k]), real(axis_phi.weight[l]), real(pars['core_sld']), 
     67                                    real(axis_theta.value[at]), real(axis_phi.value[p]), real(thickness.value[th]), 
     68                                    real(length.value[l]), real(radius.value[r]), real(pars['scale']), 
     69                                    real(radius.weight[r]), real(length.weight[l]), real(thickness.weight[th]), 
     70                                    real(axis_theta.weight[at]), real(axis_phi.weight[p]), real(pars['core_sld']), 
    7171                                    real(pars['shell_sld']), real(pars['solvent_sld']),np.uint32(size), 
    7272                                    np.uint32(self.qx.size)) 
    … …  
    7474 
    7575                            sum += self.res 
    76                             vol += radius.weight[i]*length.weight[j]*thickness.weight[f]*pow(radius.value[i]+thickness.value[f],2)\ 
    77                                    *(length.value[j]+2.0*thickness.value[f]) 
    78                             norm_vol += radius.weight[i]*length.weight[j]*thickness.weight[k] 
    79                             norm += radius.weight[i]*length.weight[j]*thickness.weight[f]*axis_theta.weight[k]\ 
    80                                     *axis_phi.weight[l] 
     76                            vol += radius.weight[r]*length.weight[l]*thickness.weight[th]*pow(radius.value[r]+thickness.value[th],2)\ 
     77                                   *(length.value[l]+2.0*thickness.value[th]) 
     78                            norm_vol += radius.weight[r]*length.weight[l]*thickness.weight[th] 
     79                            norm += radius.weight[r]*length.weight[l]*thickness.weight[th]*axis_theta.weight[at]\ 
     80                                    *axis_phi.weight[p] 
    8181 
    8282        if size>1: 

  • code_cylinder.py

    r8a21ba3 r09e15be  
    7979                        norm += radius.weight[i]*length.weight[j]*cyl_theta.weight[k]*cyl_phi.weight[l] 
    8080 
    81         if size > 1: 
    82             norm /= math.asin(1.0) 
     81       # if size > 1: 
     82        #    norm /= math.asin(1.0) 
    8383        if vol != 0.0 and norm_vol != 0.0: 
    8484            sum *= norm_vol/vol 

  • code_ellipse.py

    r8a21ba3 r09e15be  
    8484        # factor to account for the sin(theta) term in the 
    8585        # integration (see documentation). 
    86         if size > 1: 
    87             norm /= math.asin(1.0) 
    88         if vol.any() != 0.0 and norm_vol.any() != 0.0: 
     86 
     87    #    if size > 1: 
     88     #       norm /= math.asin(1.0) 
     89        if vol != 0.0 and norm_vol != 0.0: 
    8990            sum *= norm_vol/vol 
    9091 

  • code_triaxialellipse.py

    r496b252 r09e15be  
    7878                                norm += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]*theta.weight[t]*phi.weight[i]*psi.weight[s] 
    7979 
    80         if size > 1: 
    81             norm /= asin(1.0) 
     80      #  if size > 1: 
     81       #     norm /= asin(1.0) 
    8282 
    8383        if vol != 0.0 and norm_vol != 0.0: 

  • compare.py

    r8a21ba3 r09e15be  
    7777    plt.show() 
    7878 
    79 def ellipse(N=1): 
     79def ellipse(N=4): 
    8080    import sys 
    8181    import matplotlib.pyplot as plt 
    … …  
    8383    if len(sys.argv) > 1: 
    8484        N = int(sys.argv[1]) 
    85     data = load_data('JUN03289.DAT') 
     85    data = load_data('DEC07133.DAT') 
    8686    set_beam_stop(data, 0.004) 
    8787 
    … …  
    114114    plt.show() 
    115115 
     116def coreshell(N=4): 
     117    import sys 
     118    import matplotlib.pyplot as plt 
     119 
     120    if len(sys.argv) > 1: 
     121        N = int(sys.argv[1]) 
     122    data = load_data('DEC07133.DAT') 
     123    set_beam_stop(data, 0.004) 
     124 
     125    pars = dict(scale= 1.77881e-06, radius=325, thickness=25, length=34.2709, 
     126                 core_sld=1e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 
     127                 background=223.827, axis_theta=90, axis_phi=0, 
     128                 axis_theta_pd=15.8, 
     129                 radius_pd=0.1, radius_pd_n=1, radius_pd_nsigma=0, 
     130                 length_pd=0.1, length_pd_n=1, length_pd_nsigma=0, 
     131                 thickness_pd=0.1, thickness_pd_n=1, thickness_pd_nsigma=0, 
     132                 axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
     133                 axis_phi_pd=0.0008748, axis_phi_pd_n=10, axis_phi_pd_nsigma=3,) 
     134 
     135    model = sasview_model('CoreShellCylinder', **pars) 
     136    tic() 
     137    for i in range(N): 
     138        cpu = sasview_eval(model, data) 
     139    cpu_time = toc()*1000./N 
     140 
     141    from code_coreshellcyl import GpuCoreShellCylinder 
     142    model = SasModel(data, GpuCoreShellCylinder, dtype='f', **pars) 
     143    tic() 
     144    for i in range(N): 
     145        gpu = model.theory() 
     146    gpu_time = toc()*1000./N 
     147 
     148    relerr = (gpu - cpu)/cpu 
     149    print "max(|(ocl-omp)/ocl|)", max(abs(relerr)) 
     150    print "omp t=%.1f ms"%cpu_time 
     151    print "ocl t=%.1f ms"%gpu_time 
     152 
     153    plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time) 
     154    plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time) 
     155    plt.subplot(133); plot_data(data, 1e8*relerr); plt.title("relerr x 10^8"); plt.colorbar() 
     156    plt.show() 
    116157 
    117158if __name__ == "__main__": 
    118     ellipse() 
     159    coreshell() 
    119160 
    120161 

  • fit.py

    r8a21ba3 r09e15be  
    99from code_capcyl import GpuCapCylinder 
    1010from code_triaxialellipse import GpuTriEllipse 
    11 from sasmodel import SasModel, load_data, set_beam_stop 
     11from sasmodel import SasModel, load_data, set_beam_stop, set_half 
    1212 
    1313 
    14 data = load_data('JUN03289.DAT') 
    15 set_beam_stop(data, 0.004) 
     14data = load_data('DEC07133.DAT') 
     15set_beam_stop(data, 0.0048)#, outer=0.025) 
     16#set_half(data, 'right') 
    1617""" 
    17 model = SasModel(data, GpuCylinder, scale=1, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=0, 
    18                               cyl_theta=0, cyl_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,length_pd=0.1, 
    19                               length_pd_n=5, length_pd_nsigma=3, cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
    20                               cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3, dtype='float') 
    21 model.radius.range(0,100) 
     18model = SasModel(data, GpuCylinder, scale=0.09402, radius=106, length=442.298, 
     19                                    sldCyl=.291e-6, sldSolv=5.77e-6, background=0, 
     20                                    cyl_theta=0, cyl_phi=0, 
     21                                    radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=0, 
     22                                    length_pd=0.1, length_pd_n=5, length_pd_nsigma=0, 
     23                                    cyl_theta_pd=0.1, cyl_theta_pd_n=10, cyl_theta_pd_nsigma=3, 
     24                                    cyl_phi_pd=0.1, cyl_phi_pd_n=4, cyl_phi_pd_nsigma=3, 
     25                                    dtype='float') 
     26model.radius.range(0, 1000) 
    2227model.length.range(0, 1000) 
    23 model.cyl_theta.range(0,90) 
    24 model.cyl_phi.range(0,90) 
     28#model.cyl_theta.range(0,90) 
     29#model.cyl_phi.range(0,90) 
     30model.scale.range(0, 1) 
    2531""" 
     32""" 
     33model = SasModel(data, GpuEllipse, 
     34                 scale=0.0007, radius_a=42.8251, radius_b=700, sldEll=.291e-6, sldSolv=7.105e-6, 
     35                 background=10, axis_theta=0, axis_phi=0, 
     36                 radius_a_pd=0.222296, radius_a_pd_n=1, radius_a_pd_nsigma=0, 
     37                 radius_b_pd=.000128, radius_b_pd_n=1, radius_b_pd_nsigma=0, 
     38                 axis_theta_pd=24.8059, axis_theta_pd_n=40, axis_theta_pd_nsigma=3, 
     39                 axis_phi_pd=2.63698e-05, axis_phi_pd_n=20, axis_phi_pd_nsigma=3, 
     40                 dtype='float') 
    2641 
     42model.radius_a.range(15, 1000) 
     43#model.radius_b.range(15, 1000) 
     44#model.axis_theta_pd.range(0, 360) 
     45#model.background.range(0,1000) 
     46model.scale.range(0, 1) 
     47""" 
     48""" 
     49model = SasModel(data, GpuLamellar, scale=0.0131, bi_thick=41, sld_bi=.291e-6, sld_sol=5.77e-6, background=18.6, 
     50                 bi_thick_pd= 0.000009, bi_thick_pd_n=35, bi_thick_pd_nsigma=3, dtype='float') 
     51#model.bi_thick.range(0, 1000) 
     52model.scale.range(0, 1) 
     53#model.bi_thick_pd.range(0, 1000) 
     54model.background.range(0, 1000) 
     55""" 
     56""" 
     57model = SasModel(data, GpuCoreShellCylinder, 
     58                 scale= 1.77881e-06, radius=325, thickness=25, length=34.2709, 
     59                 core_sld=1e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 
     60                 background=223.827, axis_theta=90, axis_phi=0, 
     61                 axis_theta_pd=15.8, 
    2762 
    28 model = SasModel(data, GpuEllipse, scale=.027, radius_a=60, radius_b=180, sldEll=.297e-6, sldSolv=5.773e-6, background=4.9, 
    29                  axis_theta=0, axis_phi=90, radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3, radius_b_pd=0.1, radius_b_pd_n=10, 
    30                  radius_b_pd_nsigma=3, axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3, axis_phi_pd=0.1, 
    31                  axis_phi_pd_n=6, axis_phi_pd_nsigma=3, dtype='float') 
     63                 radius_pd=0.1, radius_pd_n=1, radius_pd_nsigma=0, 
     64                 length_pd=0.1, length_pd_n=1, length_pd_nsigma=0, 
     65                 thickness_pd=0.1, thickness_pd_n=1, thickness_pd_nsigma=0, 
     66                 axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
     67                 axis_phi_pd=0.0008748, axis_phi_pd_n=10, axis_phi_pd_nsigma=3, 
     68                 dtype='float') 
    3269 
    33  
    34 """ 
    35 model = SasModel(data, GpuLamellar, scale=1, bi_thick=100, sld_bi=.291e-6, sld_sol=5.77e-6, background=0, 
    36                  bi_thick_pd=0.1, bi_thick_pd_n=35, bi_thick_pd_nsigma=3, dtype='float') 
    37 """ 
    38  
    39 """ 
    40 model = SasModel(data, GpuCoreShellCylinder, scale=1, radius=64.1, thickness=1, length=266.96, core_sld=1e-6, shell_sld=4e-6, 
    41                  solvent_sld=1e-6, background=0, axis_theta=0, axis_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
    42                  length_pd=0.1, length_pd_n=10, length_pd_nsigma=3, thickness_pd=0.1, thickness_pd_n=2, thickness_pd_nsigma=3, 
    43                  axis_theta_pd=0.1, axis_theta_pd_n=2, axis_theta_pd_nsigma=3, axis_phi_pd=0.1, axis_phi_pd_n=2, 
    44                  axis_phi_pd_nsigma=3, dtype='float') 
     70model.radius.range(15, 1000) 
     71#model.length.range(0, 1000) 
     72#model.thickness.range(20, 50) 
     73#model.axis_phi.range(0, 90) 
     74#model.radius_pd.range(0, 1) 
     75#model.radius_b_pd.range(0, 1) 
     76#model.axis_theta_pd.range(0, 360) 
     77#model.axis_phi_pd.range(0, 360) 
     78#model.background.range(0,1000) 
     79model.scale.range(0, 1) 
    4580""" 
    4681 
    … …  
    5388""" 
    5489 
    55 """ 
    56 model = SasModel(data, GpuTriEllipse, scale=1, axisA=35, axisB=100, axisC=400, sldEll=1e-6, sldSolv=6.3e-6, 
    57                  background=0, theta=57, phi=57, psi=0, theta_pd=.1, theta_pd_n=4, 
    58                  theta_pd_nsigma=3, phi_pd=.1, phi_pd_n=4, phi_pd_nsigma=3, psi_pd=.1, psi_pd_n=4, psi_pd_nsigma=3, 
    59                  axisA_pd=.1, axisA_pd_n=4, axisA_pd_nsigma=3, axisB_pd=.1, axisB_pd_n=4, axisB_pd_nsigma=3, 
    60                  axisC_pd=.1, axisC_pd_n=4, axisC_pd_nsigma=3, dtype='float') 
    61 """ 
     90 
     91model = SasModel(data, GpuTriEllipse, 
     92                 scale=0.00621, axisA=20, axisB=300, axisC=900, 
     93                 sldEll=7.105e-6, sldSolv=.291e-6, 
     94                 background=15, theta=90, phi=0, psi=0, 
     95                 theta_pd=30, theta_pd_n=40, theta_pd_nsigma=3, 
     96 
     97                 phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
     98                 psi_pd=30, psi_pd_n=1, psi_pd_nsigma=0, 
     99                 axisA_pd=.1, axisA_pd_n=1, axisA_pd_nsigma=0, 
     100                 axisB_pd=.1, axisB_pd_n=1, axisB_pd_nsigma=0, 
     101                 axisC_pd=.1, axisC_pd_n=1, axisC_pd_nsigma=0, dtype='float') 
     102#model.axisA.range(15, 1000) 
     103#model.axisB.range(15, 1000) 
     104#model.axisC.range(15, 1000) 
     105#model.background.range(0,1000) 
     106model.scale.range(0, 1) 
     107model.theta_pd.range(0, 360) 
     108#model.phi_pd.range(0, 360) 
     109#model.psi_pd.range(0, 360) 
    62110 
    63111 
    64 model.scale.range(0,1) 
     112 
     113 
    65114 
    66115problem = FitProblem(model) 

  • load.py

    r5378e40 r09e15be  
    4343def demo(): 
    4444 
    45     data = load_data('JUN03289.DAT') 
     45    data = load_data('NOV07090.DAT') 
    4646    """ 
    4747    print data 

  • sasmodel.py

    r8faffcd r09e15be  
    22# -*- coding: utf-8 -*- 
    33 
     4import sys 
    45import numpy as np 
    56import pyopencl as cl 
    67from bumps.names import Parameter 
    78from sans.dataloader.loader import Loader 
    8 from sans.dataloader.manipulations import Ringcut 
     9from sans.dataloader.manipulations import Ringcut, Boxcut 
    910 
    1011 
    … …  
    1718 
    1819 
    19 def set_beam_stop(data, radius): 
     20def set_beam_stop(data, radius, outer=None): 
    2021    data.mask = Ringcut(0, radius)(data) 
     22    if outer is not None: 
     23        data.mask += Ringcut(outer,np.inf)(data) 
     24 
     25def set_half(data, half): 
     26    if half == 'left': 
     27        data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data) 
     28    if half == 'right': 
     29        data.mask += Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data) 
    2130 
    2231 
    23 def plot_data(data, iq): 
     32 
     33def plot_data(data, iq, vmin=None, vmax=None): 
    2434    from numpy.ma import masked_array 
    2535    import matplotlib.pyplot as plt 
    … …  
    2939    plt.imshow(img.reshape(128,128), 
    3040               interpolation='nearest', aspect=1, origin='upper', 
    31                extent=[xmin, xmax, ymin, ymax]) 
     41               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
    3242 
    3343 
    3444def plot_result(data, theory): 
    3545    import matplotlib.pyplot as plt 
    36     plt.subplot(1,3,1) 
    37     plot_data(data, data.data) 
    38     plt.subplot(1,3,2) 
    39     plot_data(data, theory) 
    40     plt.subplot(1,3,3) 
     46    plt.subplot(1, 3, 1) 
     47    #print "not a number",sum(np.isnan(data.data)) 
     48    #data.data[data.data<0.05] = 0.5 
     49    logdata = np.log10(data.data) 
     50    #print data.data.min(), data.data.max() 
     51    clean = logdata[~np.isnan(logdata)] 
     52    vmin, vmax = clean.min(), clean.max() 
     53    vmin, vmax = np.percentile(clean, 5), 1.05*vmax 
     54    #vmin,vmax = None,None 
     55    plot_data(data, logdata, vmin=vmin, vmax=vmax) 
     56    plt.colorbar() 
     57    plt.subplot(1, 3, 2) 
     58    plot_data(data, np.log10(theory), vmin=vmin, vmax=vmax) 
     59    plt.colorbar() 
     60    plt.subplot(1, 3, 3) 
    4161    plot_data(data, (theory-data.data)/data.err_data) 
    4262    plt.colorbar() 
    … …  
    6383    def __init__(self, data, model, dtype='float32', **kw): 
    6484        self.__dict__['_parameters'] = {} 
    65         self.index = data.mask==0 
     85        self.index = (data.mask==0) & (~np.isnan(data.data)) 
    6686        self.iq = data.data[self.index] 
    6787        self.diq = data.err_data[self.index] 
Note: See TracChangeset for help on using the changeset viewer.