Changeset a42fec0 in sasmodels


Ignore:
Timestamp:
Aug 4, 2014 5:20:07 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:
8cdb9f1
Parents:
099e053
Message:

Speed-up of 3X, compare.py working

Files:
18 edited

Legend:

Unmodified
Added
Removed

  • Kernel/Kernel-CapCyl.cpp

    r099e053 ra42fec0  
    4242        answer/=sin(alpha); 
    4343 
    44         _ptvalue[i] = rad_cyl_weight*length_weight*rad_cap_weight*theta_weight*phi_weight*vol_i*answer; 
     44        _ptvalue[i] += rad_cyl_weight*length_weight*rad_cap_weight*theta_weight*phi_weight*vol_i*answer; 
    4545     //   if (size>1) { 
    4646   //         _ptvalue[i] *= fabs(cos(thet*pi/180.0)); 

  • Kernel/Kernel-CoreShellCylinder.cpp

    r099e053 ra42fec0  
    4444        answer *= answer/(pi*(radius+thickness)*(radius+thickness)*(length+2.0*thickness))*1.0e8*scale; 
    4545 
    46         _ptvalue[i] = radius_weight*length_weight*thickness_weight*theta_weight*phi_weight*answer*pow(radius+thickness,2)*(length+2.0*thickness); 
     46        _ptvalue[i] += radius_weight*length_weight*thickness_weight*theta_weight*phi_weight*answer*pow(radius+thickness,2)*(length+2.0*thickness); 
    4747     //   if (size>1) { 
    4848       // _ptvalue[i] *= fabs(cos(axis_theta*pi/180.0)); 

  • Kernel/Kernel-Cylinder.cpp

    r099e053 ra42fec0  
    4242        real answer = sub*sub*form*acos(-1.0)*rr*rr*h*1.0e8*scale; 
    4343 
    44         _ptvalue[i] = radius_weight*length_weight*theta_weight*phi_weight*answer*pow(rr,2)*h; 
     44        _ptvalue[i] += radius_weight*length_weight*theta_weight*phi_weight*answer*pow(rr,2)*h; 
     45 
    4546       // if (size>1) { 
    4647         //   _ptvalue[i] *= fabs(cos(cyl_theta*pi/180.0)); 

  • Kernel/Kernel-Ellipse.cpp

    r099e053 ra42fec0  
    2020         ret*=ret*9.0*sub*sub*4.0/3.0*acos(-1.0)*radius_b*radius_b*radius_a*scale*(1.0e8); 
    2121 
    22          _ptvalue[i] = radius_a_weight*radius_b_weight*axis_theta_weight*radius_a*axis_phi_weight*ret*pow(radius_b, 2); 
     22         _ptvalue[i] += radius_a_weight*radius_b_weight*axis_theta_weight*radius_a*axis_phi_weight*ret*pow(radius_b, 2); 
    2323         //if(size > 1){ 
    2424          //  _ptvalue[i] *= fabs(cos(axis_theta*pi/180.0)); 

  • Kernel/Kernel-TriaxialEllipse.cpp

    r099e053 ra42fec0  
    3232        answer*=answer*sub*sub*4.0*pi/3.0*axisA*axisB*axisC*1.0e8*scale; 
    3333 
    34         _ptvalue[i] = axisA_weight*axisB_weight*axisC_weight*theta_weight*phi_weight*psi_weight*answer*axisA*axisB*axisC; 
     34        _ptvalue[i] += axisA_weight*axisB_weight*axisC_weight*theta_weight*phi_weight*psi_weight*answer*axisA*axisB*axisC; 
    3535       // if (size>1) 
    3636         //   _ptvalue[i] *= fabs(cos(theta*pi/180.0)); 

  • Models/code_capcyl.py

    rca6c007 ra42fec0  
    33 
    44import numpy as np 
    5 from math import asin, sqrt, fabs, atan 
     5from math import sqrt, fabs, atan 
    66import pyopencl as cl 
    77 
    88from weights import GaussianDispersion 
    9 from sasmodel import card 
     9from sasmodel import card, set_precision 
    1010 
    1111 
    12 def set_precision(src, qx, qy, dtype): 
    13     qx = np.ascontiguousarray(qx, dtype=dtype) 
    14     qy = np.ascontiguousarray(qy, dtype=dtype) 
    15     if np.dtype(dtype) == np.dtype('float32'): 
    16         header = """\ 
    17 #define real float 
    18 """ 
    19     else: 
    20         header = """\ 
    21 #pragma OPENCL EXTENSION cl_khr_fp64: enable 
    22 #define real double 
    23 """ 
    24     return header+src, qx, qy 
    2512 
    2613class GpuCapCylinder(object): 
    … …  
    5037 
    5138        _ctx,queue = card() 
     39        self.res[:] = 0 
     40        cl.enqueue_copy(queue, self.res_b, self.res) 
    5241        rad_cyl,length,rad_cap,theta,phi = \ 
    5342            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) 
    … …  
    7968                                        real(rad_cyl.weight[i]), real(length.weight[j]), real(theta.weight[k]), np.uint32(self.qx.size), np.uint32(size)) 
    8069 
    81                             cl.enqueue_copy(queue, self.res, self.res_b) 
    82  
    83                             sum += self.res 
    8470                            vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*vol_i 
    8571                            norm_vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m] 
    8672                            norm += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*theta.weight[k]*phi.weight[l] 
    8773 
    88         if size > 1: 
    89             norm /= asin(1.0) 
    90  
     74        #if size > 1: 
     75         #   norm /= asin(1.0) 
     76        cl.enqueue_copy(queue, self.res, self.res_b) 
     77        sum += self.res 
    9178        if vol != 0.0 and norm_vol != 0.0: 
    9279            sum *= norm_vol/vol 

  • Models/code_coreshellcyl.py

    rca6c007 ra42fec0  
    66 
    77from weights import GaussianDispersion 
    8 from sasmodel import card 
    9  
    10  
    11 def set_precision(src, qx, qy, dtype): 
    12     qx = np.ascontiguousarray(qx, dtype=dtype) 
    13     qy = np.ascontiguousarray(qy, dtype=dtype) 
    14     if np.dtype(dtype) == np.dtype('float32'): 
    15         header = """\ 
    16 #define real float 
    17 """ 
    18     else: 
    19         header = """\ 
    20 #pragma OPENCL EXTENSION cl_khr_fp64: enable 
    21 #define real double 
    22 """ 
    23     return header+src, qx, qy 
     8from sasmodel import card, set_precision 
    249 
    2510class GpuCoreShellCylinder(object): 
    … …  
    4631 
    4732        _ctx,queue = card() 
     33        self.res[:] = 0 
     34        cl.enqueue_copy(queue, self.res_b, self.res) 
    4835        radius, length, thickness, axis_phi, axis_theta = [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) 
    4936                                     for base in GpuCoreShellCylinder.PD_PARS] 
    … …  
    5643 
    5744        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0 
    58  
    59         print radius.value 
    60         print thickness.weight 
    61         print axis_phi.weight 
    62         print axis_theta.weight 
    63         print length.value 
    64  
    6545        size = len(axis_theta.weight) 
    6646 
    … …  
    7959                                    real(pars['shell_sld']), real(pars['solvent_sld']),np.uint32(size), 
    8060                                    np.uint32(self.qx.size)) 
    81                             cl.enqueue_copy(queue, self.res, self.res_b) 
    82  
    83                             sum += self.res 
    8461                            vol += radius.weight[r]*length.weight[l]*thickness.weight[th]*pow(radius.value[r]+thickness.value[th],2)\ 
    8562                                   *(length.value[l]+2.0*thickness.value[th]) 
    … …  
    9067        #if size>1: 
    9168         #   norm /= math.asin(1.0) 
     69        cl.enqueue_copy(queue, self.res, self.res_b) 
     70        sum = self.res 
    9271        if vol != 0.0 and norm_vol != 0.0: 
    9372            sum *= norm_vol/vol 

  • Models/code_cylinder.py

    rca6c007 ra42fec0  
    66 
    77from weights import GaussianDispersion 
    8 from sasmodel import card 
     8from sasmodel import card, set_precision, set_precision_1d 
    99 
    10  
    11 def set_precision(src, qx, qy, dtype): 
    12     qx = np.ascontiguousarray(qx, dtype=dtype) 
    13     qy = np.ascontiguousarray(qy, dtype=dtype) 
    14     if np.dtype(dtype) == np.dtype('float32'): 
    15         header = """\ 
    16 #define real float 
    17 """ 
    18     else: 
    19         header = """\ 
    20 #pragma OPENCL EXTENSION cl_khr_fp64: enable 
    21 #define real double 
    22 """ 
    23     return header+src, qx, qy 
    24  
    25 def set_precision_1d(src, q, dtype): 
    26     q = np.ascontiguousarray(q, dtype=dtype) 
    27     if np.dtype(dtype) == np.dtype('float32'): 
    28         header = """\ 
    29 #define real float 
    30 """ 
    31     else: 
    32         header = """\ 
    33 #pragma OPENCL EXTENSION cl_khr_fp64: enable 
    34 #define real double 
    35 """ 
    36     return header+src, q 
    3710 
    3811class GpuCylinder(object): 
    … …  
    5528        self.qx_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx) 
    5629        self.qy_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy) 
    57         self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, qx.nbytes) 
     30        self.res_b = cl.Buffer(ctx, cl.mem_flags.READ_WRITE, self.qx.nbytes) 
    5831        self.res = np.empty_like(self.qx) 
    5932 
    … …  
    6134 
    6235        _ctx,queue = card() 
     36        self.res[:] = 0 
     37        cl.enqueue_copy(queue, self.res_b, self.res) 
    6338        radius, length, cyl_theta, cyl_phi = \ 
    6439            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) 
    … …  
    6843        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True) 
    6944        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True) 
    70         cyl_theta.value, cyl_theta.weight = cyl_theta.get_weights(pars['cyl_theta'], -90, 180, False) 
    71         cyl_phi.value, cyl_phi.weight = cyl_phi.get_weights(pars['cyl_phi'], -90, 180, False) 
     45        cyl_theta.value, cyl_theta.weight = cyl_theta.get_weights(pars['cyl_theta'], -np.inf, np.inf, False) 
     46        cyl_phi.value, cyl_phi.weight = cyl_phi.get_weights(pars['cyl_phi'], -np.inf, np.inf, False) 
    7247 
    7348        #Perform the computation, with all weight points 
    … …  
    8762                                           real(cyl_phi.weight[l]), real(cyl_theta.value[k]), real(cyl_phi.value[l]), 
    8863                                           np.uint32(self.qx.size), np.uint32(size)) 
    89                         cl.enqueue_copy(queue, self.res, self.res_b) 
    90                         sum += self.res 
     64 
    9165                        vol += radius.weight[i]*length.weight[j]*pow(radius.value[i], 2)*length.value[j] 
    9266                        norm_vol += radius.weight[i]*length.weight[j] 
    … …  
    9569       # if size > 1: 
    9670        #    norm /= math.asin(1.0) 
     71        cl.enqueue_copy(queue, self.res, self.res_b) 
     72        sum = self.res 
    9773        if vol != 0.0 and norm_vol != 0.0: 
    9874            sum *= norm_vol/vol 

  • Models/code_ellipse.py

    rca6c007 ra42fec0  
    66 
    77from weights import GaussianDispersion 
    8 from sasmodel import card 
    9  
    10  
    11 def set_precision(src, qx, qy, dtype): 
    12     qx = np.ascontiguousarray(qx, dtype=dtype) 
    13     qy = np.ascontiguousarray(qy, dtype=dtype) 
    14     if np.dtype(dtype) == np.dtype('float32'): 
    15         header = """\ 
    16 #define real float 
    17 """ 
    18     else: 
    19         header = """\ 
    20 #pragma OPENCL EXTENSION cl_khr_fp64: enable 
    21 #define real double 
    22 """ 
    23     return header+src, qx, qy 
     8from sasmodel import card, set_precision 
    249 
    2510class GpuEllipse(object): 
    … …  
    4631    #b_n = radius_b # want, a_n = radius_a # want, etc 
    4732        _ctx,queue = card() 
     33        self.res[:] = 0 
     34        cl.enqueue_copy(queue, self.res_b, self.res) 
    4835        radius_a, radius_b, axis_theta, axis_phi = \ 
    4936            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) 
    … …  
    7663                                        real(axis_phi.value[l]), self.qx_b, self.qy_b, self.res_b, 
    7764                                        np.uint32(self.qx.size), np.uint32(len(axis_theta.weight))) 
    78                         #copy result back from buffer 
    79                         cl.enqueue_copy(queue, self.res, self.res_b) 
    80                         sum += self.res 
     65 
    8166                        vol += radius_a.weight[i]*radius_b.weight[j]*pow(radius_b.value[j], 2)*radius_a.value[i] 
    8267                        norm_vol += radius_a.weight[i]*radius_b.weight[j] 
    … …  
    8873    #    if size > 1: 
    8974     #       norm /= math.asin(1.0) 
     75        cl.enqueue_copy(queue, self.res, self.res_b) 
     76        sum += self.res 
    9077        if vol != 0.0 and norm_vol != 0.0: 
    9178            sum *= norm_vol/vol 

  • Models/code_lamellar.py

    rca6c007 ra42fec0  
    55import pyopencl as cl 
    66from Models.weights import GaussianDispersion 
    7  
    8 def set_precision(src, qx, qy, dtype): 
    9     qx = np.ascontiguousarray(qx, dtype=dtype) 
    10     qy = np.ascontiguousarray(qy, dtype=dtype) 
    11     if dtype == 'double': 
    12         header = """\ 
    13 #define real float 
    14 """ 
    15     else: 
    16         header = """\ 
    17 #pragma OPENCL EXTENSION cl_khr_fp64: enable 
    18 #define real double 
    19 """ 
    20     return header+src, qx, qy 
     7from sasmodel import set_precision 
    218 
    229 

  • Models/code_triaxialellipse.py

    rca6c007 ra42fec0  
    66 
    77from weights import GaussianDispersion 
    8 from sasmodel import card 
    9  
    10  
    11 def set_precision(src, qx, qy, dtype): 
    12     qx = np.ascontiguousarray(qx, dtype=dtype) 
    13     qy = np.ascontiguousarray(qy, dtype=dtype) 
    14     if np.dtype(dtype) == np.dtype('float32'): 
    15         header = """\ 
    16 #define real float 
    17 """ 
    18     else: 
    19         header = """\ 
    20 #pragma OPENCL EXTENSION cl_khr_fp64: enable 
    21 #define real double 
    22 """ 
    23     return header+src, qx, qy 
     8from sasmodel import card, set_precision 
    249 
    2510class GpuTriEllipse: 
    … …  
    4530 
    4631        _ctx,queue = card() 
     32        self.res[:] = 0 
     33        cl.enqueue_copy(queue, self.res_b, self.res) 
    4734        axisA, axisB, axisC, theta, phi, psi = \ 
    4835            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) 
    … …  
    7259                                            real(axisA.weight[a]), real(axisB.weight[b]), real(axisC.weight[c]), real(psi.weight[s]), 
    7360                                            real(phi.weight[i]), real(theta.weight[t]), np.uint32(self.qx.size), np.uint32(size)) 
    74                                 cl.enqueue_copy(queue, self.res, self.res_b) 
    75                                 sum += self.res 
    7661 
    7762                                vol += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]*axisA.value[a]*axisB.value[b]*axisC.value[c] 
    … …  
    8166      #  if size > 1: 
    8267       #     norm /= asin(1.0) 
    83  
     68        cl.enqueue_copy(queue, self.res, self.res_b) 
     69        sum = self.res 
    8470        if vol != 0.0 and norm_vol != 0.0: 
    8571            sum *= norm_vol/vol 

  • Models/weights.py

    rdbb0048 ra42fec0  
    1919        sigma = width * center if relative else width 
    2020        if sigma == 0: 
    21             return np.array([center, 1.], 'd') 
     21            return np.array([center],'d'), np.array([1.], 'd') 
    2222        x = center + np.linspace(-nsigmas * sigma, +nsigmas * sigma, npts) 
    2323        x = x[(x >= min) & (x <= max)] 

  • Testers/test.py

    r099e053 ra42fec0  
     1impor 
     2 
     3 
     4 
     5 
    16import numpy as np 
    27 
    … …  
    1217 
    1318print sum 
     19 

  • compare.py

    r473183c ra42fec0  
    44import datetime 
    55 
    6 from Models.sasmodel import SasModel, load_data, set_beam_stop, plot_data 
     6from sasmodel import SasModel, load_data, set_beam_stop, plot_data 
    77 
    88 
    … …  
    4646    if len(sys.argv) > 1: 
    4747        N = int(sys.argv[1]) 
    48     data = load_data('JUN03289.DAT') 
     48    data = load_data('December/DEC07098.DAT') 
    4949    set_beam_stop(data, 0.004) 
    5050 
     51    dtype = "float" 
    5152    pars = dict( 
    52         scale=1, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=0, 
    53         cyl_theta=0, cyl_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,length_pd=0.1, 
    54         length_pd_n=5, length_pd_nsigma=3, cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
     53        scale=.08, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1, 
     54        cyl_theta=0, cyl_phi=0, 
     55        radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
     56        length_pd=0.1,length_pd_n=5, length_pd_nsigma=3, 
     57        cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
    5558        cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3, 
    5659        ) 
    … …  
    6265    cpu_time = toc()*1000./N 
    6366 
    64     from code_cylinder import GpuCylinder 
     67 
     68    from Models.code_cylinder import GpuCylinder 
    6569    model = SasModel(data, GpuCylinder, dtype='f', **pars) 
    6670    tic() 
    … …  
    6973    gpu_time = toc()*1000./N 
    7074 
     75    print "gpu/cpu", max(gpu/cpu) 
     76    cpu *= max(gpu/cpu) 
    7177    relerr = (gpu - cpu)/cpu 
    7278    print "max(|(ocl-omp)/ocl|)", max(abs(relerr)) 
    … …  
    7682    plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time) 
    7783    plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time) 
    78     plt.subplot(133); plot_data(data, 1e8*relerr); plt.title("relerr x 10^8"); plt.colorbar() 
     84    plt.subplot(133); plot_data(data, relerr); plt.title("relerr x 10^8"); plt.colorbar() 
    7985    plt.show() 
    8086 
    … …  
    8591    if len(sys.argv) > 1: 
    8692        N = int(sys.argv[1]) 
    87     data = load_data('DEC07133.DAT') 
     93    data = load_data('DEC07106.DAT') 
    8894    set_beam_stop(data, 0.004) 
    8995 
    … …  
    99105    cpu_time = toc()*1000./N 
    100106 
    101     from code_ellipse import GpuEllipse 
     107    from Models.code_ellipse import GpuEllipse 
    102108    model = SasModel(data, GpuEllipse, dtype='f', **pars) 
    103109    tic() 
    … …  
    116122    plt.show() 
    117123 
    118 def coreshell(N=4): 
     124def coreshell(N=1): 
    119125    import sys 
    120126    import matplotlib.pyplot as plt 
    … …  
    122128    if len(sys.argv) > 1: 
    123129        N = int(sys.argv[1]) 
    124     data = load_data('DEC07133.DAT') 
     130    data = load_data('December/DEC07098.DAT') 
    125131    set_beam_stop(data, 0.004) 
    126132 
    … …  
    141147    cpu_time = toc()*1000./N 
    142148 
    143     from code_coreshellcyl import GpuCoreShellCylinder 
     149    from Models.code_coreshellcyl import GpuCoreShellCylinder 
    144150    model = SasModel(data, GpuCoreShellCylinder, dtype='f', **pars) 
    145151    tic() 

  • fit.py

    r099e053 ra42fec0  
    33 
    44from bumps.names import * 
    5 from sasmodel import SasModel, load_data, set_beam_stop 
     5from sasmodel import SasModel, load_data, set_beam_stop, set_half, set_top 
    66from Models.code_capcyl import GpuCapCylinder 
    77from Models.code_coreshellcyl import GpuCoreShellCylinder 
    … …  
    1313""" IMPORT THE DATA USED """ 
    1414#data = load_data('December/Tangential/Sector0/DEC07133.ABS') 
    15 data = load_data('December/DEC07235.DAT') 
     15data = load_data('December/DEC07102.DAT') 
    1616 
    1717""" SET INNER BEAM STOP, OUTER RING, AND MASK HALF OF THE DATA """ 
    18 set_beam_stop(data, 0.0052)#, outer=0.025) 
    19 #set_half(data, 'left') 
     18set_beam_stop(data, 0.00669)#, outer=0.025) 
     19set_top(data, -.018) 
     20#set_half(data, 'right') 
    2021 
    2122""" 
    … …  
    2627length=1000, 
    2728background=21, 
    28 sldCyl=.291e-6,sldSolv=5.77e-6, 
     29sldCyl=.291e-6,sldSolv=7.105e-6, 
    2930radius_pd=0.1,radius_pd_n=10,radius_pd_nsigma=0, 
    3031length_pd=0.1,length_pd_n=5,length_pd_nsigma=0, 
    … …  
    3536model = SasModel(data, GpuEllipse, 
    3637scale=0.08, 
    37 radius_a=15, radius_b=200, 
     38radius_a=15, radius_b=800, 
    3839sldEll=.291e-6, sldSolv=7.105e-6, 
    3940background=0, 
    40 axis_theta=0, axis_phi=0, 
    41 axis_theta_pd=20, axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
     41axis_theta=90, axis_phi=0, 
     42axis_theta_pd=15, axis_theta_pd_n=40, axis_theta_pd_nsigma=3, 
    4243 
    4344radius_a_pd=0.222296, radius_a_pd_n=1, radius_a_pd_nsigma=0, 
    … …  
    4849 
    4950# SET THE FITTING PARAMETERS 
    50 #model.radius_a.range(15, 1000) 
    51 #model.radius_b.range(15, 1000) 
    52 model.axis_theta_pd.range(0, 360) 
     51model.radius_a.range(15, 1000) 
     52model.radius_b.range(15, 1000) 
     53#model.axis_theta_pd.range(0, 360) 
    5354#model.background.range(0,1000) 
    5455#model.scale.range(0, 1) 
    5556""" 
     57 
    5658""" 
    57  
    5859model = SasModel(data, GpuLamellar, 
    59 scale=0.70, 
    60 bi_thick=5, 
    61 sld_bi=.291e-6,sld_sol=5.77e-6, 
    62 background=85.23, 
    63 bi_thick_pd= 0.0013, bi_thick_pd_n=5, bi_thick_pd_nsigma=3, 
     60scale=0.08, 
     61bi_thick=19.2946, 
     62sld_bi=5.38e-6,sld_sol=7.105e-6, 
     63background=0.003, 
     64bi_thick_pd= 0.37765, bi_thick_pd_n=10, bi_thick_pd_nsigma=3, 
    6465dtype='float') 
    6566 
    6667# SET THE FITTING PARAMETERS 
    67 model.bi_thick.range(0, 1000) 
    68 model.scale.range(0, 1) 
     68#model.bi_thick.range(0, 1000) 
     69#model.scale.range(0, 1) 
    6970#model.bi_thick_pd.range(0, 1000) 
    7071#model.background.range(0, 1000) 
     72model.sld_bi.range(0, 1) 
     73 
    7174""" 
    72  
    73 model = SasModel(data, GpuCylinder, 
    74 scale=0.08, radius=46, length=719, 
    75 sldCyl=.291e-6, sldSolv=5.77e-6, background=0, 
    76 cyl_theta=0, cyl_phi=0, 
    77 cyl_theta_pd=23, cyl_theta_pd_n=10, cyl_theta_pd_nsigma=3, 
     75if 1: 
     76    model = SasModel(data, GpuCylinder, 
     77    scale=0.0104, 
     78    radius=92.5, 
     79    length=798.3, 
     80    sldCyl=.29e-6, 
     81    sldSolv=7.105e-6, 
     82    background=5, 
     83    cyl_theta=0, 
     84    cyl_phi=0, 
     85    cyl_theta_pd=22.11, 
     86    cyl_theta_pd_n=20, 
     87    cyl_theta_pd_nsigma=3, 
     88    radius_pd=.0084, 
     89    radius_pd_n=10, 
     90    radius_pd_nsigma=3, 
     91    length_pd=0.493, 
     92    length_pd_n=10, 
     93    length_pd_nsigma=3, 
     94    cyl_phi_pd=0, 
     95    cyl_phi_pd_n=1, 
     96    cyl_phi_pd_nsigma=3, 
     97    dtype='float') 
    7898 
    7999 
    80 radius_pd=0.1, radius_pd_n=20, radius_pd_nsigma=0, 
    81 length_pd=0.1, length_pd_n=20, length_pd_nsigma=0, 
    82 cyl_phi_pd=0.1, cyl_phi_pd_n=20, cyl_phi_pd_nsigma=0, 
    83 dtype='float') 
    84100 
    85 # SET THE FITTING PARAMETERS 
    86 model.radius.range(0, 25) 
    87 #model.length.range(0, 1000) 
    88 #model.cyl_theta_pd.range(0,90) 
    89 #model.scale.range(0, 1) 
    90 #model.background.range(0, 1000) 
     101    # SET THE FITTING PARAMETERS 
     102    #model.radius.range(1, 500) 
     103    #model.length.range(1, 4000) 
     104    #model.cyl_theta.range(-90,100) 
     105    #model.cyl_theta_pd.range(0, 90) 
     106    #model.cyl_theta_pd_n = model.cyl_theta_pd + 5 
     107    #model.radius_pd.range(0, 90) 
     108    #model.length_pd.range(0, 90) 
     109    model.scale.range(0, 1) 
     110    #model.background.range(0, 100) 
     111    #model.sldCyl.range(0, 1) 
    91112 
    92113 
    93 """ 
    94 model = SasModel(data, GpuCoreShellCylinder, 
    95                  scale= 0.08, radius=200, thickness=30, length=2000, 
    96                  core_sld=7e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 
    97                  background=0, axis_theta=0, axis_phi=0, 
     114if 0: 
     115    model = SasModel(data, GpuCoreShellCylinder, 
     116                     scale= .00031, radius=19.5, thickness=30, length=22, 
     117                     core_sld=7.105e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 
     118                     background=0.2, axis_theta=0, axis_phi=0, 
    98119 
    99                  radius_pd=0.38, radius_pd_n=10, radius_pd_nsigma=3, 
    100                  length_pd=.9, length_pd_n=10, length_pd_nsigma=3, 
    101                  thickness_pd=0.1, thickness_pd_n=1, thickness_pd_nsigma=0, 
    102                  axis_theta_pd=10, axis_theta_pd_n=40, axis_theta_pd_nsigma=3, 
    103                  axis_phi_pd=0.1, axis_phi_pd_n=1, axis_phi_pd_nsigma=0, 
    104                  dtype='float') 
     120                     radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3, 
     121                     length_pd=0.26, length_pd_n=10, length_pd_nsigma=3, 
     122                     thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=1, 
     123                     axis_theta_pd=1, axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
     124                     axis_phi_pd=0.1, axis_phi_pd_n=1, axis_phi_pd_nsigma=1, 
     125                     dtype='float') 
    105126 
    106 # SET THE FITTING PARAMETERS 
    107 model.radius.range(15, 1000) 
    108 #model.length.range(0, 1000) 
    109 #model.thickness.range(20, 50) 
    110 #model.axis_phi.range(0, 90) 
    111 #model.radius_pd.range(0, 1) 
    112 #model.radius_b_pd.range(0, 1) 
    113 #model.axis_theta_pd.range(0, 360) 
    114 #model.axis_phi_pd.range(0, 360) 
    115 #model.background.range(0,1000) 
    116 model.scale.range(0, 1) 
    117 """ 
     127    # SET THE FITTING PARAMETERS 
     128    #model.radius.range(115, 1000) 
     129    #model.length.range(0, 2500) 
     130    #model.thickness.range(18, 38) 
     131    #model.thickness_pd.range(0, 1) 
     132    #model.axis_phi.range(0, 90) 
     133    #model.radius_pd.range(0, 1) 
     134    #model.length_pd.range(0, 1) 
     135    #model.axis_theta_pd.range(0, 360) 
     136    #model.background.range(0,5) 
     137    model.scale.range(0, 1) 
     138 
    118139 
    119140 
    … …  
    135156""" 
    136157model = SasModel(data, GpuTriEllipse, 
    137                  scale=0.0036, axisA=118, axisB=70, axisC=800, 
     158                 scale=0.08, axisA=15, axisB=20, axisC=500, 
    138159                 sldEll=7.105e-6, sldSolv=.291e-6, 
    139                  background=15, theta=90, phi=0, psi=0, 
    140                  theta_pd=22, theta_pd_n=40, theta_pd_nsigma=3, 
     160                 background=5, theta=0, phi=0, psi=0, 
     161                 theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3, 
    141162                 phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    142163                 psi_pd=30, psi_pd_n=1, psi_pd_nsigma=0, 
    … …  
    150171#model.axisC.range(15, 1000) 
    151172#model.background.range(0,1000) 
    152 model.scale.range(0, 1) 
    153173#model.theta_pd.range(0, 360) 
    154174#model.phi_pd.range(0, 360) 

  • fit2.py

    r473183c ra42fec0  
    33 
    44from bumps.names import * 
    5 from cylcode import GpuCylinder 
    6 from Models.sasmodel import SasModel, load_data, set_beam_stop 
     5from Models.code_coreshellcyl import GpuCoreShellCylinder 
     6from sasmodel import SasModel, load_data, set_beam_stop, set_top 
    77 
    88 
    9 radial_data = load_data('JUN03289.DAT') 
    10 set_beam_stop(radial_data, 0.004) 
    11 trans_data = load_data('JUN03305.DAT') 
    12 set_beam_stop(trans_data, 0.004) 
     9radial_data = load_data('December/DEC07267.DAT') 
     10set_beam_stop(radial_data, 0.00669, outer=0.025) 
     11set_top(radial_data, -.0185) 
     12 
     13tan_data = load_data('December/DEC07266.DAT') 
     14set_beam_stop(tan_data, 0.00669, outer=0.025) 
     15set_top(tan_data, -.0185) 
    1316 
    1417 
    … …  
    1720dtype='float32' 
    1821radial = SasModel(radial_data, 
    19                   GpuCylinder, dtype=dtype, 
    20                   scale=1, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=0, 
    21                   cyl_theta=0, cyl_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,length_pd=0.1, 
    22                   length_pd_n=5, length_pd_nsigma=3, cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
    23                   cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3) 
    24 trans = SasModel(trans_data, 
    25                   GpuCylinder, dtype=dtype, 
     22                 GpuCoreShellCylinder, 
     23                 scale= 3.75e-7, radius=378, thickness=30, length=1806, 
     24                 core_sld=7.105e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 
     25                 background=0.2, axis_theta=0, axis_phi=90, 
     26 
     27                 radius_pd=0.26, radius_pd_n=20, radius_pd_nsigma=3, 
     28                 length_pd=0.26, length_pd_n=20, length_pd_nsigma=3, 
     29                 thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=0, 
     30                 axis_theta_pd=1, axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
     31                 axis_phi_pd=0.1, axis_phi_pd_n=1, axis_phi_pd_nsigma=0, 
     32                 dtype='float') 
     33tan = SasModel(tan_data, 
     34                  GpuCoreShellCylinder, dtype=dtype, 
    2635                  **radial.parameters()) 
    2736 
    28 radial.radius.range(0,100) 
    29 radial.length.range(0, 1000) 
    30 radial.cyl_theta.range(0,90) 
    31 radial.cyl_phi.range(0,90) 
    32 radial.scale.range(0,10) 
    33 trans.cyl_theta = radial.cyl_theta + 90. 
     37radial.radius.range(15, 1000) 
     38radial.length.range(0, 2500) 
     39#radial.thickness.range(18, 38) 
     40#radial.thickness_pd.range(0, 1) 
     41#radial.axis_phi.range(0, 90) 
     42#radial.radius_pd.range(0, 1) 
     43#radial.length_pd.range(0, 1) 
     44#radial.axis_theta_pd.range(0, 360) 
     45#radial.background.range(0,5) 
     46#radial.scale.range(0, 1) 
     47radial.axis_phi = tan.axis_phi + 90 
    3448 
    3549 
    36 problem = FitProblem([radial,trans]) 
     50problem = FitProblem([radial,tan]) 
    3751 
    3852 

  • multi-fit.py

    rd772f5d ra42fec0  
    33 
    44from bumps.names import * 
    5 from code_lamellar import GpuLamellar 
    6 from code_ellipse import GpuEllipse 
     5from Models.code_lamellar import GpuLamellar 
     6from Models.code_ellipse import GpuEllipse 
     7from Models.code_cylinder import GpuCylinder 
    78from multisasmodels import SasModel, load_data, set_beam_stop, set_half 
    89import numpy as np 
    910 
    10 data = load_data('DEC07282.DAT') 
     11data = load_data('December/DEC07238.DAT') 
    1112set_beam_stop(data, 0.0052)#, outer=0.025) 
    1213#set_half(data, 'left') 
    1314 
    14 truth = SasModel(data, GpuEllipse, 
    15                 scale=.0011, radius_a=45.265, radius_b=600.8, sldEll=.291e-6, sldSolv=7.105e-6, 
    16                  background=8.30161, axis_theta=0, axis_phi=0, 
    17                  radius_a_pd=0.222296, radius_a_pd_n=1, radius_a_pd_nsigma=0, 
    18                  radius_b_pd=.000128, radius_b_pd_n=1, radius_b_pd_nsigma=0, 
    19                  axis_theta_pd=20, axis_theta_pd_n=40, axis_theta_pd_nsigma=3, 
    20                  axis_phi_pd=2.63698e-05, axis_phi_pd_n=20, axis_phi_pd_nsigma=0, 
    21                  dtype='float') 
     15truth = SasModel(data, GpuCylinder, 
     16scale=0.08, radius=46, length=719, 
     17sldCyl=.291e-6, sldSolv=5.77e-6, background=0, 
     18cyl_theta=90, cyl_phi=0, 
     19cyl_theta_pd=23, cyl_theta_pd_n=40, cyl_theta_pd_nsigma=3, 
     20radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
     21length_pd=0.1, length_pd_n=10, length_pd_nsigma=3, 
     22cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3, 
     23dtype='float') 
     24 
    2225 
    2326#model.radius_a.range(15, 1000) 
    … …  
    3033 
    3134 
    32 lies = SasModel(data, GpuLamellar, scale=0, bi_thick=5, sld_bi=.291e-6, sld_sol=5.77e-6, background=85.23, 
    33                  bi_thick_pd= 0.0013, bi_thick_pd_n=5, bi_thick_pd_nsigma=3, dtype='float') 
     35lies = SasModel(data, GpuCylinder, 
     36scale=0.08, radius=46, length=719, 
     37sldCyl=.291e-6, sldSolv=5.77e-6, background=0, 
     38cyl_theta=90, cyl_phi=0, 
     39cyl_theta_pd=23, cyl_theta_pd_n=40, cyl_theta_pd_nsigma=3, 
     40radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
     41length_pd=0.1, length_pd_n=10, length_pd_nsigma=3, 
     42cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3, 
     43dtype='float') 
     44 
    3445 
    3546 
    … …  
    4354 
    4455 
    45 a = np.add(np.multiply(arrayone, .5), np.multiply(arraytwo, 0)) 
     56a = np.add(np.multiply(arrayone, .08), np.multiply(arraytwo, .08)) 
    4657truth.set_result(a) 
    4758 

  • sasmodel.py

    r099e053 ra42fec0  
    1616    return data 
    1717 
     18def set_precision(src, qx, qy, dtype): 
     19    qx = np.ascontiguousarray(qx, dtype=dtype) 
     20    qy = np.ascontiguousarray(qy, dtype=dtype) 
     21    if np.dtype(dtype) == np.dtype('float32'): 
     22        header = """\ 
     23#define real float 
     24""" 
     25    else: 
     26        header = """\ 
     27#pragma OPENCL EXTENSION cl_khr_fp64: enable 
     28#define real double 
     29""" 
     30    return header+src, qx, qy 
     31 
     32def set_precision_1d(src, q, dtype): 
     33    q = np.ascontiguousarray(q, dtype=dtype) 
     34    if np.dtype(dtype) == np.dtype('float32'): 
     35        header = """\ 
     36#define real float 
     37""" 
     38    else: 
     39        header = """\ 
     40#pragma OPENCL EXTENSION cl_khr_fp64: enable 
     41#define real double 
     42""" 
     43    return header+src, q 
    1844 
    1945def set_beam_stop(data, radius, outer=None): 
    … …  
    2854 
    2955def set_half(data, half): 
     56    if half == 'right': 
     57        data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data) 
    3058    if half == 'left': 
    31         data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data) 
    32     if half == 'right': 
    3359        data.mask += Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data) 
    3460 
    35  
     61def set_top(data, max): 
     62    data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data) 
    3663 
    3764def plot_data(data, iq, vmin=None, vmax=None): 
    … …  
    6996    plt.colorbar() 
    7097    plt.subplot(1, 3, 3) 
     98    print abs(mresid).max() 
    7199    plot_data(data, mresid) 
    72100    plt.colorbar() 
    … …  
    110138    def __init__(self, data, model, dtype='float32', **kw): 
    111139        self.__dict__['_parameters'] = {} 
     140        #self.name = data.filename 
    112141        self.is2D = hasattr(data,'qx_data') 
    113142        self.data = data 
    … …  
    135164        pars.update((base+'_pd_nsigma', 3) for base in model.PD_PARS) 
    136165        pars.update(kw) 
    137         self._parameters = dict((k, Parameter.default(v, name=k)) for k, v in pars.items()) 
     166        for k,v in pars.items(): 
     167            setattr(self, k, Parameter.default(v, name=k)) 
     168        self._parameter_names = set(pars.keys()) 
     169        self.update() 
     170 
     171    def update(self): 
     172        self._cache = {} 
    138173 
    139174    def numpoints(self): 
    … …  
    141176 
    142177    def parameters(self): 
    143         return self._parameters 
    144  
    145     def __getattr__(self, par): 
    146         return self._parameters[par] 
    147  
    148     def __setattr__(self, par, val): 
    149         if par in self._parameters: 
    150             self._parameters[par] = val 
    151         else: 
    152             self.__dict__[par] = val 
     178        return dict((k,getattr(self,k)) for k in self._parameter_names) 
    153179 
    154180    def theory(self): 
    155         pars = dict((k,v.value) for k,v in self._parameters.items()) 
    156         result = self.gpu.eval(pars) 
    157         return result 
     181        if 'theory' not in self._cache: 
     182            pars = dict((k,getattr(self,k).value) for k in self._parameter_names) 
     183            #print pars 
     184            self._cache['theory'] = self.gpu.eval(pars) 
     185        return self._cache['theory'] 
    158186 
    159187    def residuals(self): 
    … …  
    178206        pass 
    179207 
    180     def update(self): 
    181         pass 
    182  
    183  
    184208def demo(): 
    185209    data = load_data('DEC07086.DAT') 
Note: See TracChangeset for help on using the changeset viewer.