Changeset 1726b21 in sasmodels


Ignore:
Timestamp:
Aug 8, 2014 1:45:51 PM (11 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:
ae7d639
Parents:
8cdb9f1
Message:

cylinder now MUCH faster!

Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

  • Models/code_capcyl.py

    ra42fec0 r1726b21  
    1212 
    1313class GpuCapCylinder(object): 
    14     PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'length':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0, 
     14    PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'len_cyl':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0, 
    1515             'theta':0, 'phi':0} 
    1616 
    17     PD_PARS = ['rad_cyl', 'length', 'rad_cap', 'theta', 'phi'] 
     17    PD_PARS = ['rad_cyl', 'len_cyl', 'rad_cap', 'theta', 'phi'] 
    1818 
    1919    def __init__(self, qx, qy, dtype='float32'): 
    … …  
    3939        self.res[:] = 0 
    4040        cl.enqueue_copy(queue, self.res_b, self.res) 
    41         rad_cyl,length,rad_cap,theta,phi = \ 
     41 
     42        rad_cyl,len_cyl,rad_cap,theta,phi = \ 
    4243            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) 
    4344             for base in GpuCapCylinder.PD_PARS] 
    … …  
    4546        rad_cyl.value, rad_cyl.weight = rad_cyl.get_weights(pars['rad_cyl'], 0, 10000, True) 
    4647        rad_cap.value, rad_cap.weight = rad_cap.get_weights(pars['rad_cap'], 0, 10000, True) 
    47         length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True) 
     48        len_cyl.value, len_cyl.weight = len_cyl.get_weights(pars['len_cyl'], 0, 10000, True) 
    4849        theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False) 
    4950        phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False) 
    … …  
    5657        for i in xrange(len(rad_cyl.weight)): 
    5758            for m in xrange(len(rad_cap.weight)): 
    58                 for j in xrange(len(length.weight)): 
     59                for j in xrange(len(len_cyl.weight)): 
     60 
     61                    hDist = -1.0*sqrt(fabs(rad_cap.value[m]*rad_cap.value[m]-rad_cyl.value[i]*rad_cyl.value[i])) 
     62                    vol_i = 4.0*atan(1.0)*rad_cyl.value[i]*rad_cyl.value[i]*len_cyl.value[j]+2.0*4.0*atan(1.0)/3.0\ 
     63                                    *((rad_cap.value[m]-hDist)*(rad_cap.value[m]-hDist)*(2*rad_cap.value[m]+hDist)) 
     64                    vol += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]*vol_i 
     65                    norm_vol += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m] 
     66 
    5967                    for k in xrange(len(theta.weight)): 
    6068                        for l in xrange(len(phi.weight)): 
    61                             hDist = -1.0*sqrt(fabs(rad_cap.value[m]*rad_cap.value[m]-rad_cyl.value[i]*rad_cyl.value[i])) 
    62                             vol_i = 4.0*atan(1.0)*rad_cyl.value[i]*rad_cyl.value[i]*length.value[j]+2.0*4.0*atan(1.0)/3.0\ 
    63                                             *((rad_cap.value[m]-hDist)*(rad_cap.value[m]-hDist)*(2*rad_cap.value[m]+hDist)) 
     69 
    6470                            self.prg.CapCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, 
    65                                         real(vol_i), real(hDist), real(rad_cyl.value[i]), real(rad_cap.value[m]), real(length.value[j]), 
     71                                        real(vol_i), real(hDist), real(rad_cyl.value[i]), real(rad_cap.value[m]), real(len_cyl.value[j]), 
    6672                                        real(theta.value[k]), real(phi.value[l]), real(sub), real(pars['scale']), 
    6773                                        real(phi.weight[l]), real(theta.weight[k]), real(rad_cap.weight[m]), 
    68                                         real(rad_cyl.weight[i]), real(length.weight[j]), real(theta.weight[k]), np.uint32(self.qx.size), np.uint32(size)) 
     74                                        real(rad_cyl.weight[i]), real(len_cyl.weight[j]), real(theta.weight[k]), np.uint32(self.qx.size), np.uint32(size)) 
    6975 
    70                             vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*vol_i 
    71                             norm_vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m] 
    72                             norm += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*theta.weight[k]*phi.weight[l] 
     76                            norm += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]*theta.weight[k]*phi.weight[l] 
     77 
     78 
    7379 
    7480        #if size > 1: 

  • Models/code_coreshellcyl.py

    ra42fec0 r1726b21  
    4848        for r in xrange(len(radius.weight)): 
    4949            for l in xrange(len(length.weight)): 
    50                 for at in xrange(len(axis_theta.weight)): 
    51                     for p in xrange(len(axis_phi.weight)): 
    52                         for th in xrange(len(thickness.weight)): 
     50                for th in xrange(len(thickness.weight)): 
     51 
     52                    vol += radius.weight[r]*length.weight[l]*thickness.weight[th]*pow(radius.value[r]+thickness.value[th],2)\ 
     53                                   *(length.value[l]+2.0*thickness.value[th]) 
     54                    norm_vol += radius.weight[r]*length.weight[l]*thickness.weight[th] 
     55 
     56                    for at in xrange(len(axis_theta.weight)): 
     57                        for p in xrange(len(axis_phi.weight)): 
    5358 
    5459                            self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, 
    … …  
    5964                                    real(pars['shell_sld']), real(pars['solvent_sld']),np.uint32(size), 
    6065                                    np.uint32(self.qx.size)) 
    61                             vol += radius.weight[r]*length.weight[l]*thickness.weight[th]*pow(radius.value[r]+thickness.value[th],2)\ 
    62                                    *(length.value[l]+2.0*thickness.value[th]) 
    63                             norm_vol += radius.weight[r]*length.weight[l]*thickness.weight[th] 
     66 
    6467                            norm += radius.weight[r]*length.weight[l]*thickness.weight[th]*axis_theta.weight[at]\ 
    6568                                    *axis_phi.weight[p] 

  • Models/code_cylinder.py

    ra42fec0 r1726b21  
    66 
    77from weights import GaussianDispersion 
    8 from sasmodel import card, set_precision, set_precision_1d 
     8from sasmodel import card, set_precision, set_precision_1d, tic, toc 
    99 
    1010 
    … …  
    3333    def eval(self, pars): 
    3434 
     35        tic() 
    3536        _ctx,queue = card() 
    3637        self.res[:] = 0 
    … …  
    5556        for i in xrange(len(radius.weight)): 
    5657            for j in xrange(len(length.weight)): 
     58 
     59                vol += radius.weight[i]*length.weight[j]*pow(radius.value[i], 2)*length.value[j] 
     60                norm_vol += radius.weight[i]*length.weight[j] 
     61 
    5762                for k in xrange(len(cyl_theta.weight)): 
    5863                    for l in xrange(len(cyl_phi.weight)): 
    … …  
    6368                                           np.uint32(self.qx.size), np.uint32(size)) 
    6469 
    65                         vol += radius.weight[i]*length.weight[j]*pow(radius.value[i], 2)*length.value[j] 
    66                         norm_vol += radius.weight[i]*length.weight[j] 
    6770                        norm += radius.weight[i]*length.weight[j]*cyl_theta.weight[k]*cyl_phi.weight[l] 
     71 
    6872 
    6973       # if size > 1: 
    … …  
    7478            sum *= norm_vol/vol 
    7579 
     80        print toc()*1000, self.qx.shape[0] 
    7681        return sum/norm+pars['background'] 
    7782 

  • Models/code_ellipse.py

    ra42fec0 r1726b21  
    5353            for j in xrange(len(radius_b.weight)): 
    5454                #Average over theta distribution 
     55                vol += radius_a.weight[i]*radius_b.weight[j]*pow(radius_b.value[j], 2)*radius_a.value[i] 
     56                norm_vol += radius_a.weight[i]*radius_b.weight[j] 
     57 
    5558                for k in xrange(len(axis_theta.weight)): 
    5659                    #Average over phi distribution 
    … …  
    6467                                        np.uint32(self.qx.size), np.uint32(len(axis_theta.weight))) 
    6568 
    66                         vol += radius_a.weight[i]*radius_b.weight[j]*pow(radius_b.value[j], 2)*radius_a.value[i] 
    67                         norm_vol += radius_a.weight[i]*radius_b.weight[j] 
    6869                        norm += radius_a.weight[i]*radius_b.weight[j]*axis_theta.weight[k]*axis_phi.weight[l] 
     70 
    6971        # Averaging in theta needs an extra normalization 
    7072        # factor to account for the sin(theta) term in the 

  • Models/code_triaxialellipse.py

    ra42fec0 r1726b21  
    99 
    1010class GpuTriEllipse: 
    11     PARS = {'scale':1, 'axisA':35, 'axisB':100, 'axisC':400, 'sldEll':1e-6, 'sldSolv':6.3e-6, 'background':0, 
    12             'theta':0, 'phi':0, 'psi':0} 
     11    PARS = {'scale':1, 'semi_axisA':35, 'semi_axisB':100, 'semi_axisC':400, 'sldEll':1e-6, 'sldSolv':6.3e-6, 'background':0, 
     12            'axis_theta':0, 'axis_phi':0, 'axis_psi':0} 
    1313 
    14     PD_PARS = ['axisA', 'axisB', 'axisC', 'theta', 'phi', 'psi'] 
     14    PD_PARS = ['semi_axisA', 'semi_axisB', 'semi_axisC', 'axis_theta', 'axis_phi', 'axis_psi'] 
    1515 
    1616    def __init__(self, qx, qy, dtype='float32'): 
    … …  
    3232        self.res[:] = 0 
    3333        cl.enqueue_copy(queue, self.res_b, self.res) 
    34         axisA, axisB, axisC, theta, phi, psi = \ 
     34        semi_axisA, semi_axisB, semi_axisC, axis_theta, axis_phi, axis_psi = \ 
    3535            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) 
    3636             for base in GpuTriEllipse.PD_PARS] 
    3737 
    38         axisA.value, axisA.weight = axisA.get_weights(pars['axisA'], 0, 10000, True) 
    39         axisB.value, axisB.weight = axisB.get_weights(pars['axisB'], 0, 10000, True) 
    40         axisC.value, axisC.weight = axisC.get_weights(pars['axisC'], 0, 10000, True) 
    41         theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False) 
    42         phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False) 
    43         psi.value, psi.weight = psi.get_weights(pars['psi'], -90, 180, False) 
     38        semi_axisA.value, semi_axisA.weight = semi_axisA.get_weights(pars['semi_axisA'], 0, 10000, True) 
     39        semi_axisB.value, semi_axisB.weight = semi_axisB.get_weights(pars['semi_axisB'], 0, 10000, True) 
     40        semi_axisC.value, semi_axisC.weight = semi_axisC.get_weights(pars['semi_axisC'], 0, 10000, True) 
     41        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False) 
     42        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False) 
     43        axis_psi.value, axis_psi.weight = axis_psi.get_weights(pars['axis_psi'], -90, 180, False) 
    4444 
    4545        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0 
    46         size = len(theta.weight) 
     46        size = len(axis_theta.weight) 
    4747        sub = pars['sldEll'] - pars['sldSolv'] 
    4848 
    4949        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64 
    50         for a in xrange(len(axisA.weight)): 
    51             for b in xrange(len(axisB.weight)): 
    52                 for c in xrange(len(axisC.weight)): 
    53                     for t in xrange(len(theta.weight)): 
    54                         for i in xrange(len(phi.weight)): 
    55                             for s in xrange(len(psi.weight)): 
     50        for a in xrange(len(semi_axisA.weight)): 
     51            for b in xrange(len(semi_axisB.weight)): 
     52                for c in xrange(len(semi_axisC.weight)): 
     53 
     54                    vol += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]*semi_axisA.value[a]*semi_axisB.value[b]*semi_axisC.value[c] 
     55                    norm_vol += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c] 
     56 
     57                    for t in xrange(len(axis_theta.weight)): 
     58                        for i in xrange(len(axis_phi.weight)): 
     59                            for s in xrange(len(axis_psi.weight)): 
    5660                                self.prg.TriaxialEllipseKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, 
    57                                             real(sub), real(pars['scale']), real(axisA.value[a]), real(axisB.value[b]), 
    58                                             real(axisC.value[c]), real(phi.value[i]), real(theta.value[t]), real(psi.value[s]), 
    59                                             real(axisA.weight[a]), real(axisB.weight[b]), real(axisC.weight[c]), real(psi.weight[s]), 
    60                                             real(phi.weight[i]), real(theta.weight[t]), np.uint32(self.qx.size), np.uint32(size)) 
     61                                            real(sub), real(pars['scale']), real(semi_axisA.value[a]), real(semi_axisB.value[b]), 
     62                                            real(semi_axisC.value[c]), real(axis_phi.value[i]), real(axis_theta.value[t]), real(axis_psi.value[s]), 
     63                                            real(semi_axisA.weight[a]), real(semi_axisB.weight[b]), real(semi_axisC.weight[c]), real(axis_psi.weight[s]), 
     64                                            real(axis_phi.weight[i]), real(axis_theta.weight[t]), np.uint32(self.qx.size), np.uint32(size)) 
    6165 
    62                                 vol += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]*axisA.value[a]*axisB.value[b]*axisC.value[c] 
    63                                 norm_vol += axisA.weight[a]*axisB.weight[b]*axisC.weight[c] 
    64                                 norm += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]*theta.weight[t]*phi.weight[i]*psi.weight[s] 
     66                                norm += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]*axis_theta.weight[t]*axis_phi.weight[i]*axis_psi.weight[s] 
     67 
    6568 
    6669      #  if size > 1: 

  • compare.py

    ra42fec0 r1726b21  
    22# -*- coding: utf-8 -*- 
    33 
    4 import datetime 
    5  
    6 from sasmodel import SasModel, load_data, set_beam_stop, plot_data 
    7  
    8  
    9 TIC = None 
    10 def tic(): 
    11     global TIC 
    12     TIC = datetime.datetime.now() 
    13  
    14 def toc(): 
    15     now = datetime.datetime.now() 
    16     return (now-TIC).total_seconds() 
     4from sasmodel import SasModel, load_data, set_beam_stop, plot_data, tic, toc 
     5import numpy as np 
     6 
    177 
    188def sasview_model(modelname, **pars): 
    … …  
    4030    return theory 
    4131 
    42 def cyl(N=1): 
     32def plot(data, cpumodel, gpumodel, N, pars): 
     33 
     34    model = SasModel(data, gpumodel, dtype='f', **pars) 
     35    tic() 
     36    for i in range(N): 
     37        #pars['scale'] = np.random.rand() 
     38        model.update() 
     39        gpu = model.theory() 
     40    gpu_time = toc()*1000./N 
     41    print "ocl t=%.1f ms"%gpu_time 
     42    tic() 
     43    for i in range(N): 
     44        cpu = sasview_eval(cpumodel, data) 
     45    cpu_time = toc()*1000./N 
     46    print "omp t=%.1f ms"%cpu_time 
     47 
     48    import matplotlib.pyplot as plt 
     49 
     50    print "gpu/cpu", max(gpu/cpu) 
     51    cpu *= max(gpu/cpu) 
     52    relerr = (gpu - cpu)/cpu 
     53    print "max(|(ocl-omp)/ocl|)", max(abs(relerr)) 
     54 
     55 
     56    plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time) 
     57    plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time) 
     58    plt.subplot(133); plot_data(data, relerr); plt.title("relerr x 10^8"); plt.colorbar() 
     59    plt.show() 
     60 
     61def newlen(N): 
    4362    import sys 
    44     import matplotlib.pyplot as plt 
    4563 
    4664    if len(sys.argv) > 1: 
    … …  
    4967    set_beam_stop(data, 0.004) 
    5068 
     69    return data, N 
     70 
     71def cyl(N=1): 
     72 
    5173    dtype = "float" 
    5274    pars = dict( 
    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, 
     75        scale=.003, radius=64.1, length=66.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1, 
     76        cyl_theta=90, cyl_phi=0, 
    5577        radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
    5678        length_pd=0.1,length_pd_n=5, length_pd_nsigma=3, 
    5779        cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
    58         cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3, 
    59         ) 
    60  
     80        cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3) 
     81 
     82    from Models.code_cylinder_f import GpuCylinder as gpumodel 
    6183    model = sasview_model('Cylinder', **pars) 
    62     tic() 
    63     for i in range(N): 
    64         cpu = sasview_eval(model, data) 
    65     cpu_time = toc()*1000./N 
    66  
    67  
    68     from Models.code_cylinder import GpuCylinder 
    69     model = SasModel(data, GpuCylinder, dtype='f', **pars) 
    70     tic() 
    71     for i in range(N): 
    72         gpu = model.theory() 
    73     gpu_time = toc()*1000./N 
    74  
    75     print "gpu/cpu", max(gpu/cpu) 
    76     cpu *= max(gpu/cpu) 
    77     relerr = (gpu - cpu)/cpu 
    78     print "max(|(ocl-omp)/ocl|)", max(abs(relerr)) 
    79     print "omp t=%.1f ms"%cpu_time 
    80     print "ocl t=%.1f ms"%gpu_time 
    81  
    82     plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time) 
    83     plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time) 
    84     plt.subplot(133); plot_data(data, relerr); plt.title("relerr x 10^8"); plt.colorbar() 
    85     plt.show() 
    86  
    87 def ellipse(N=4): 
    88     import sys 
    89     import matplotlib.pyplot as plt 
    90  
    91     if len(sys.argv) > 1: 
    92         N = int(sys.argv[1]) 
    93     data = load_data('DEC07106.DAT') 
    94     set_beam_stop(data, 0.004) 
     84    data, N = newlen(N) 
     85 
     86    plot(data, model, gpumodel, N, pars) 
     87 
     88 
     89def ellipse(N=1): 
    9590 
    9691    pars = dict(scale=.027, radius_a=60, radius_b=180, sldEll=.297e-6, sldSolv=5.773e-6, background=4.9, 
    97                 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, 
    98                 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, 
    99                 axis_phi_pd_n=6, axis_phi_pd_nsigma=3,) 
    100  
     92                axis_theta=0, axis_phi=90, 
     93                radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3, 
     94                radius_b_pd=0.1, radius_b_pd_n=10, radius_b_pd_nsigma=3, 
     95                axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3, 
     96                axis_phi_pd=0.1, axis_phi_pd_n=6, axis_phi_pd_nsigma=3,) 
     97 
     98    from Models.code_ellipse import GpuEllipse as gpumodel 
    10199    model = sasview_model('Ellipsoid', **pars) 
    102     tic() 
    103     for i in range(N): 
    104         cpu = sasview_eval(model, data) 
    105     cpu_time = toc()*1000./N 
    106  
    107     from Models.code_ellipse import GpuEllipse 
    108     model = SasModel(data, GpuEllipse, dtype='f', **pars) 
    109     tic() 
    110     for i in range(N): 
    111         gpu = model.theory() 
    112     gpu_time = toc()*1000./N 
    113  
    114     relerr = (gpu - cpu)/cpu 
    115     print "max(|(ocl-omp)/ocl|)", max(abs(relerr)) 
    116     print "omp t=%.1f ms"%cpu_time 
    117     print "ocl t=%.1f ms"%gpu_time 
    118  
    119     plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time) 
    120     plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time) 
    121     plt.subplot(133); plot_data(data, 1e8*relerr); plt.title("relerr x 10^8"); plt.colorbar() 
    122     plt.show() 
     100 
     101    data, N = newlen(N) 
     102    plot(data, model, gpumodel, N, pars) 
    123103 
    124104def coreshell(N=1): 
    125     import sys 
    126     import matplotlib.pyplot as plt 
    127  
    128     if len(sys.argv) > 1: 
    129         N = int(sys.argv[1]) 
    130     data = load_data('December/DEC07098.DAT') 
    131     set_beam_stop(data, 0.004) 
     105 
     106    data, N = newlen(N) 
    132107 
    133108    pars = dict(scale= 1.77881e-06, radius=325, thickness=25, length=34.2709, 
    … …  
    135110                 background=223.827, axis_theta=90, axis_phi=0, 
    136111                 axis_theta_pd=15.8, 
    137                  radius_pd=0.1, radius_pd_n=1, radius_pd_nsigma=0, 
    138                  length_pd=0.1, length_pd_n=1, length_pd_nsigma=0, 
    139                  thickness_pd=0.1, thickness_pd_n=1, thickness_pd_nsigma=0, 
    140                  axis_theta_pd_n=20, axis_theta_pd_nsigma=3, 
    141                  axis_phi_pd=0.0008748, axis_phi_pd_n=60, axis_phi_pd_nsigma=3,) 
     112                 radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
     113                 length_pd=0.1, length_pd_n=10, length_pd_nsigma=3, 
     114                 thickness_pd=0.1, thickness_pd_n=5, thickness_pd_nsigma=3, 
     115                 axis_theta_pd_n=20, axis_theta_pd_nsigma=5, 
     116                 axis_phi_pd=0.0008748, axis_phi_pd_n=5, axis_phi_pd_nsigma=3,) 
    142117 
    143118    model = sasview_model('CoreShellCylinder', **pars) 
    144     tic() 
    145     for i in range(N): 
    146         cpu = sasview_eval(model, data) 
    147     cpu_time = toc()*1000./N 
    148  
    149     from Models.code_coreshellcyl import GpuCoreShellCylinder 
    150     model = SasModel(data, GpuCoreShellCylinder, dtype='f', **pars) 
    151     tic() 
    152     for i in range(N): 
    153         gpu = model.theory() 
    154     gpu_time = toc()*1000./N 
    155  
    156     relerr = (gpu - cpu)/cpu 
    157     print "max(|(ocl-omp)/ocl|)", max(abs(relerr)) 
    158     print "omp t=%.1f ms"%cpu_time 
    159     print "ocl t=%.1f ms"%gpu_time 
    160  
    161     plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time) 
    162     plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time) 
    163     plt.subplot(133); plot_data(data, 1e8*relerr); plt.title("relerr x 10^8"); plt.colorbar() 
    164     plt.show() 
     119    from Models.code_coreshellcyl import GpuCoreShellCylinder as gpumodel 
     120 
     121    plot(data, model, gpumodel, N, pars) 
     122 
     123def trellipse(N=1): 
     124 
     125    data, N = newlen(N) 
     126 
     127    pars = dict(scale=0.08, semi_axisA=15, semi_axisB=20, semi_axisC=500, 
     128                sldEll=7.105e-6, sldSolv=.291e-6, 
     129                background=5, axis_theta=0, axis_phi=0, axis_psi=0, 
     130                axis_theta_pd=20, axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
     131                axis_phi_pd=.1, axis_phi_pd_n=10, axis_phi_pd_nsigma=3, 
     132                axis_psi_pd=30, axis_psi_pd_n=5, axis_psi_pd_nsigma=3, 
     133                semi_axisA_pd=.1, semi_axisA_pd_n=5, semi_axisA_pd_nsigma=3, 
     134                semi_axisB_pd=.1, semi_axisB_pd_n=5, semi_axisB_pd_nsigma=3, 
     135                semi_axisC_pd=.1, semi_axisC_pd_n=5, semi_axisC_pd_nsigma=3,) 
     136 
     137    model = sasview_model('TriaxialEllipsoid', **pars) 
     138    from Models.code_triaxialellipse import GpuTriEllipse as gpumodel 
     139 
     140    plot(data, model, gpumodel, N, pars) 
     141 
     142def lamellar(N=1): 
     143 
     144    data, N = newlen(N) 
     145 
     146    pars = dict(scale=0.08, 
     147                bi_thick=19.2946, 
     148                sld_bi=5.38e-6,sld_sol=7.105e-6, 
     149                background=0.003, 
     150                bi_thick_pd= 0.37765, bi_thick_pd_n=40, bi_thick_pd_nsigma=3) 
     151 
     152    model = sasview_model('Lamellar', **pars) 
     153    from Models.code_lamellar import GpuLamellar as gpumodel 
     154 
     155    plot(data, model, gpumodel, N, pars) 
     156 
     157def cap(N=1): 
     158 
     159    data, N = newlen(N) 
     160 
     161    pars = dict(scale=.08, rad_cyl=20, rad_cap=40, len_cyl=400, 
     162                sld_capcyl=1e-6, sld_solv=6.3e-6, 
     163                background=0, theta=0, phi=0, 
     164                rad_cyl_pd=.1, rad_cyl_pd_n=10, rad_cyl_pd_nsigma=3, 
     165                rad_cap_pd=.1, rad_cap_pd_n=10, rad_cap_pd_nsigma=3, 
     166                len_cyl_pd=.1, len_cyl_pd_n=3, len_cyl_pd_nsigma=3, 
     167                theta_pd=.1, theta_pd_n=3, theta_pd_nsigma=3, 
     168                phi_pd=.1, phi_pd_n=3, phi_pd_nsigma=3) 
     169 
     170 
     171    model = sasview_model('CappedCylinder', **pars) 
     172    from Models.code_capcyl import GpuCapCylinder as gpumodel 
     173 
     174    plot(data, model, gpumodel, N, pars) 
     175 
    165176 
    166177if __name__ == "__main__": 
    167     coreshell() 
    168  
    169  
    170  
    171  
    172  
    173  
    174  
    175  
    176  
    177  
    178  
    179  
    180  
    181  
    182  
    183  
    184  
    185  
    186  
    187  
    188  
    189  
    190  
    191  
    192  
    193  
    194  
    195  
    196  
     178    cyl(2) 
     179 
     180 
     181 
     182 
     183 
     184 
     185 
     186 
     187 
     188 
     189 
     190 
     191 
     192 
     193 
     194 
     195 
     196 
     197 
     198 
     199 
     200 
     201 
     202 
     203 
     204 
     205 
     206 
     207 

  • fit.py

    r8cdb9f1 r1726b21  
    66from Models.code_capcyl import GpuCapCylinder 
    77from Models.code_coreshellcyl import GpuCoreShellCylinder 
    8 from Models.code_cylinder import GpuCylinder, OneDGpuCylinder 
     8from Models.code_cylinder_f import GpuCylinder 
     9#from Models.code_cylinder import GpuCylinder, OneDGpuCylinder 
    910from Models.code_ellipse import GpuEllipse 
    1011from Models.code_lamellar import GpuLamellar 
    … …  
    1314""" IMPORT THE DATA USED """ 
    1415#data = load_data('December/Tangential/Sector0/DEC07133.ABS') 
    15 data = load_data('December/DEC07102.DAT') 
     16data = load_data('December/DEC07098.DAT') 
    1617 
    1718""" SET INNER BEAM STOP, OUTER RING, AND MASK HALF OF THE DATA """ 
    18 set_beam_stop(data, 0.00669)#, outer=0.025) 
    19 set_top(data, -.018) 
     19set_beam_stop(data, 0.004)#, outer=0.025) 
     20#set_top(data, -.018) 
    2021#set_half(data, 'right') 
    2122 
    … …  
    4142    radius_b_pd=.000128, radius_b_pd_n=1, radius_b_pd_nsigma=0, 
    4243    axis_phi_pd=2.63698e-05, axis_phi_pd_n=20, axis_phi_pd_nsigma=0, 
    43     dtype='float') 
     44    dtype='float32') 
    4445 
    4546 
    … …  
    5960    background=0.003, 
    6061    bi_thick_pd= 0.37765, bi_thick_pd_n=10, bi_thick_pd_nsigma=3, 
    61     dtype='float') 
     62    dtype='float32') 
    6263 
    6364    # SET THE FITTING PARAMETERS 
    … …  
    7071 
    7172if 1: 
    72     model = SasModel(data, GpuCylinder, 
    73     scale=0.0104, radius=92.5, length=798.3, 
    74     sldCyl=.29e-6, sldSolv=7.105e-6, background=5, 
    75     cyl_theta=0, cyl_phi=0, 
    76     cyl_theta_pd=22.11, cyl_theta_pd_n=20, cyl_theta_pd_nsigma=3, 
    77     radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3, 
    78     length_pd=0.493, length_pd_n=10, length_pd_nsigma=3, 
    79     cyl_phi_pd=0, cyl_phi_pd_n=1, cyl_phi_pd_nsigma=3, 
    80     dtype='float') 
     73    """ 
     74    pars = dict(scale=0.0023, radius=92.5, length=798.3, 
     75        sldCyl=.29e-6, sldSolv=7.105e-6, background=5, 
     76        cyl_theta=0, cyl_phi=0, 
     77        cyl_theta_pd=22.11, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
     78        radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3, 
     79        length_pd=0.493, length_pd_n=10, length_pd_nsigma=3, 
     80        cyl_phi_pd=0, cyl_phi_pd_n=5, cyl_phi_pd_nsigma=3,) 
     81        """ 
     82    pars = dict( 
     83        scale=.01, radius=64.1, length=66.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1, 
     84        cyl_theta=90, cyl_phi=0, 
     85        radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
     86        length_pd=0.1,length_pd_n=5, length_pd_nsigma=3, 
     87        cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
     88        cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3) 
     89    model = SasModel(data, GpuCylinder, dtype="float32", **pars) 
    8190 
    8291 
    8392 
    8493    # SET THE FITTING PARAMETERS 
    85     #model.radius.range(1, 500) 
    86     #model.length.range(1, 4000) 
     94    model.radius.range(1, 500) 
     95    model.length.range(1, 4000) 
    8796    #model.cyl_theta.range(-90,100) 
    8897    #model.cyl_theta_pd.range(0, 90) 
    … …  
    106115    axis_theta_pd=1, axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
    107116    axis_phi_pd=0.1, axis_phi_pd_n=1, axis_phi_pd_nsigma=1, 
    108     dtype='float') 
     117    dtype='float32') 
    109118 
    110119    # SET THE FITTING PARAMETERS 
    … …  
    124133if 0: 
    125134    model = SasModel(data, GpuCapCylinder, 
    126     scale=1, rad_cyl=20, rad_cap=40, length=400, 
     135    scale=.08, rad_cyl=20, rad_cap=40, len_cyl=400, 
    127136    sld_capcyl=1e-6, sld_solv=6.3e-6, 
    128137    background=0, theta=0, phi=0, 
    129     rad_cyl_pd=.1, rad_cyl_pd_n=1, rad_cyl_pd_nsigma=0, 
    130     rad_cap_pd=.1, rad_cap_pd_n=1, rad_cap_pd_nsigma=0, 
    131     length_pd=.1, length_pd_n=1, length_pd_nsigma=0, 
     138    rad_cyl_pd=.1, rad_cyl_pd_n=5, rad_cyl_pd_nsigma=3, 
     139    rad_cap_pd=.1, rad_cap_pd_n=5, rad_cap_pd_nsigma=3, 
     140    len_cyl_pd=.1, len_cyl_pd_n=1, len_cyl_pd_nsigma=0, 
    132141    theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0, 
    133142    phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    134     dtype='float') 
     143    dtype='float32') 
    135144 
    136145    model.scale.range(0, 1) 
    … …  
    147156    axisA_pd=.1, axisA_pd_n=1, axisA_pd_nsigma=0, 
    148157    axisB_pd=.1, axisB_pd_n=1, axisB_pd_nsigma=0, 
    149     axisC_pd=.1, axisC_pd_n=1, axisC_pd_nsigma=0, dtype='float') 
     158    axisC_pd=.1, axisC_pd_n=1, axisC_pd_nsigma=0, dtype='float32') 
    150159 
    151160    # SET THE FITTING PARAMETERS 

  • sasmodel.py

    ra42fec0 r1726b21  
    22# -*- coding: utf-8 -*- 
    33 
     4import datetime 
    45import numpy as np 
    56import pyopencl as cl 
    … …  
    89from sans.dataloader.manipulations import Ringcut, Boxcut 
    910 
     11 
     12TIC = None 
     13def tic(): 
     14    global TIC 
     15    TIC = datetime.datetime.now() 
     16 
     17def toc(): 
     18    now = datetime.datetime.now() 
     19    return (now-TIC).total_seconds() 
    1020 
    1121def load_data(filename): 
    … …  
    2131    if np.dtype(dtype) == np.dtype('float32'): 
    2232        header = """\ 
     33 
     34#define REAL(x) (x##f) 
    2335#define real float 
    2436""" 
    … …  
    2638        header = """\ 
    2739#pragma OPENCL EXTENSION cl_khr_fp64: enable 
     40#define REAL(x) (x) 
    2841#define real double 
    2942""" 
Note: See TracChangeset for help on using the changeset viewer.