Changeset 8a20be5 in sasmodels


Ignore:
Timestamp:
Jul 10, 2014 1:05:08 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:
8faffcd
Parents:
5378e40
Message:

Added a fit2 (fits two different models at different angles)
(preliminary) Added CoreshellCyl? and CapCyl? Kernels
(preliminary) Updated kernels to include functions

Files:
6 added
8 edited

Legend:

Unmodified
Added
Removed

  • Kernel-CoreShellCylinder.cpp

    r5378e40 r8a20be5  
    3737        if (besarg1 == 0.0){be1 = 0.5;} 
    3838        else{ 
    39             if((ax=fabs(besarg1)) < 8.0) 
    40             { 
    41                 y=besarg1*besarg1; 
    42                 ans1=besarg1*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+y*(15704.48260+y*(-30.16036606)))))); 
    43                 ans2=144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0)))); 
    44                 ans=ans1/ans2; 
    45             } 
    46             else 
    47             { 
    48                 z=8.0/ax; 
    49                 y=z*z; 
    50                 xx=ax-2.356194491; 
    51                 ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6)))); 
    52                 ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6))); 
    53                 ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2); 
    54                 if (besarg1 < 0.0) {ans *= -1;} 
    55             } 
    56             be1 = ans/besarg1; 
     39            bel = NR_BessJ1(besarg1)/besarg1 
    5740        } 
    5841        if (besarg2 == 0.0){be2 = 0.5;} 
    … …  
    10992 
    11093 
    111  
    112  
    113  
    114  
    115  
    116  
    117  
    118  
    119  
    120  
    121  
    122  
    123  
    124  
    125  
    126  
    127  
    128  
    129  
    130  
    131  
    132  
    133  
    134  
    135  
    136  
    137  
    138  
    139  
    140  
    141  
    142  
    143  
    144  
    145  
    146  
    147  
    148  
    149  
    150  
    151  
    152  
    153  
    154  
    155  
    156  
    157  
    158  
    159  
    160  

  • Kernel-Cylinder.cpp

    r5378e40 r8a20be5  
    1 __kernel void CylinderKernel(__global const float *qx, global const float *qy, __global float *_ptvalue, const float sub, 
    2 const float rr, const float h, const float scale, const float radius_weight, const float length_weight, 
    3 const float theta_weight, const float phi_weight, const float cyl_theta, 
    4 const float cyl_phi, const int count, const int size) 
     1__kernel void CylinderKernel(__global const real *qx, global const real *qy, __global real *_ptvalue, const real sub, 
     2const real rr, const real h, const real scale, const real radius_weight, const real length_weight, 
     3const real theta_weight, const real phi_weight, const real cyl_theta, 
     4const real cyl_phi, const int count, const int size) 
    55{ 
    66        // qq is the q-value for the calculation (1/A) 
    … …  
    1111    if(i < count) 
    1212    { 
    13         float qq = sqrt(qx[i]*qx[i]+qy[i]*qy[i]); 
     13        real qq = sqrt(qx[i]*qx[i]+qy[i]*qy[i]); 
    1414 
    15         float pi = 4.0*atan(1.0); 
    16         float theta = cyl_theta*pi/180.0; 
    17         float phi = cyl_phi*pi/180.0; 
     15        real pi = 4.0*atan(1.0); 
     16        real theta = cyl_theta*pi/180.0; 
     17        real phi = cyl_phi*pi/180.0; 
    1818 
    19         float cyl_x = cos(theta)*cos(phi); 
    20         float cyl_y = sin(theta); 
    21         float cos_val = cyl_x*(qx[i]/qq) + cyl_y*(qy[i]/qq); 
     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); 
    2222 
    23         float alpha = acos(cos_val); 
     23        real alpha = acos(cos_val); 
    2424        if(alpha == 0.0){ 
    2525            alpha = 1.0e-26; 
    2626        } 
    27         float besarg = qq*rr*sin(alpha); 
    28         float siarg = qq*h/2*cos(alpha); 
     27        real besarg = qq*rr*sin(alpha); 
     28        real siarg = qq*h/2*cos(alpha); 
     29        real be=0.0; real si=0.0; 
    2930 
    30             float xx=0.0; float y=0.0; float bj=0.0; float ans1=0.0; float ans2=0.0; float z=0.0; float answer=0.0; 
    31             float contrast=0.0; float form=0.0; float be=0.0; float si=0.0; 
     31        real bj = NR_BessJ1(besarg); 
    3232 
    33         float ax = fabs(besarg); 
     33        real d1 = qq*rr*sin(alpha); 
    3434 
    35         if(ax < 8.0) { 
    36             y=besarg*besarg; 
    37             ans1=besarg*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+y*(15704.48260+y*(-30.16036606)))))); 
    38             ans2=144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0)))); 
    39             bj=ans1/ans2; 
     35        if (besarg == 0.0){ 
     36            be = sin(alpha); 
    4037        } 
    4138        else{ 
    42             z=8.0/ax; 
    43             y=z*z; 
    44             xx=ax - 2.356194491; 
    45             ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6)))); 
    46             ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6))); 
    47             bj=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2); 
    48  
    49             if(besarg < 0.0){bj*=-1;} 
     39            be = bj*bj*4.0*sin(alpha)/(d1*d1); 
     40        } 
     41        if(siarg == 0.0){ 
     42            si = 1.0; 
     43        } 
     44        else{ 
     45            si = sin(siarg)*sin(siarg)/(siarg*siarg); 
    5046        } 
    5147 
    52         float d1 = qq*rr*sin(alpha); 
    53  
    54         if (besarg == 0.0) {be = sin(alpha);} 
    55         else {be = bj*bj*4.0*sin(alpha)/(d1*d1);} 
    56         if(siarg == 0.0) {si = 1.0;} 
    57         else{si = sin(siarg)*sin(siarg)/(siarg*siarg);} 
    58  
    59         form = be*si/sin(alpha); 
    60         answer = sub*sub*form*acos(-1.0)*rr*rr*h*1.0e8*scale; 
     48        real form = be*si/sin(alpha); 
     49        real answer = sub*sub*form*acos(-1.0)*rr*rr*h*1.0e8*scale; 
    6150 
    6251        _ptvalue[i] = radius_weight*length_weight*theta_weight*phi_weight*answer*pow(rr,2)*h; 
    … …  
    6453            _ptvalue[i] *= fabs(cos(cyl_theta*pi/180.0)); 
    6554        } 
    66 } 
     55    } 
    6756} 
    6857 

  • Kernel-Lamellar.cpp

    r5378e40 r8a20be5  
    1 __kernel void LamellarKernel(__global const float *qx, global const float *qy, __global float *ret, const float bi_thick, 
    2  const float scale, const float sub, const float background, const int length) 
     1#ifndef real 
     2# define real float 
     3#endif 
     4 
     5__kernel void LamellarKernel(__global const real *qx, global const real *qy, __global real *ret, const real bi_thick, 
     6 const real scale, const real sub, const real background, const int length) 
    37{ 
    48    int i = get_global_id(0); 
    59    if(i < length) 
    610    { 
    7         float q = sqrt(qx[i]*qx[i]+qy[i]*qy[i]); 
    8         float pi = 4.0*atan(1.0); 
    9         float Pq = 2.0*sub*(sub/q)/q*(1.0-cos(q*bi_thick)); 
     11        real q = sqrt(qx[i]*qx[i]+qy[i]*qy[i]); 
     12        real pi = 4.0*atan(1.0); 
     13        real Pq = 2.0*sub*(sub/q)/q*(1.0-cos(q*bi_thick)); 
    1014        ret[i] = 2.0*pi*scale*Pq/(q*q)/bi_thick*1.0e8; 
    1115        ret[i] += background; 

  • coreshellcylcode.py

    r5378e40 r8a20be5  
    66import pyopencl as cl 
    77from weights import GaussianDispersion 
     8from sasmodel import card 
     9 
     10def set_precision(src, qx, qy, dtype): 
     11    qx = np.ascontiguousarray(qx, dtype=dtype) 
     12    qy = np.ascontiguousarray(qy, dtype=dtype) 
     13    if np.dtype(dtype) == np.dtype('float32'): 
     14        header = """\ 
     15#define real float 
     16""" 
     17    else: 
     18        header = """\ 
     19#pragma OPENCL EXTENSION cl_khr_fp64: enable 
     20#define real double 
     21""" 
     22    return header+src, qx, qy 
    823 
    924class GpuCoreShellCylinder(object): 
    10     PARS = {'scale':1, 'radius':1, 'thickness':1, 'length':1, 'core_sld':1e-6, 'shell_sld':1e-6, 'solvent_sld':0, 
     25    PARS = {'scale':1, 'radius':1, 'thickness':1, 'length':1, 'core_sld':1e-6, 'shell_sld':-1e-6, 'solvent_sld':0, 
    1126            'background':0, 'axis_theta':0, 'axis_phi':0} 
    1227    PD_PARS = ['radius', 'length', 'thickness', 'axis_phi', 'axis_theta'] 
    1328 
    14     def __init__(self, qx, qy): 
    15         self.qx = np.asarray(qx, np.float32) 
    16         self.qy = np.asarray(qy, np.float32) 
     29    def __init__(self, qx, qy, dtype='float32'): 
    1730        #create context, queue, and build program 
    18         self.ctx = cl.create_some_context() 
    19         self.queue = cl.CommandQueue(self.ctx) 
    20         self.prg = cl.Program(self.ctx, open('Kernel-CoreShellCylinder.cpp').read()).build() 
     31        ctx,_queue = card() 
     32        src, qx, qy = set_precision(open('NR_BessJ1.cpp').read()+"\n"+open('Kernel-CoreShellCylinder.cpp').read(), qx, qy, dtype=dtype) 
     33        self.prg = cl.Program(ctx, src).build() 
     34        self.qx, self.qy = qx, qy 
     35 
    2136 
    2237        #buffers 
    2338        mf = cl.mem_flags 
    24         self.qx_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx) 
    25         self.qy_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy) 
    26         self.res_b = cl.Buffer(self.ctx, mf.WRITE_ONLY, qx.nbytes) 
    27         self.res = np.empty_like(self.qx) 
     39        self.qx_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx) 
     40        self.qy_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy) 
     41        self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, qx.nbytes) 
     42        self.res = np.empty_like(qx) 
    2843 
    2944    def eval(self, pars): 
    3045 
     46        _ctx,queue = card() 
    3147        radius, length, thickness, axis_phi, axis_theta = [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) 
    3248                                     for base in GpuCoreShellCylinder.PD_PARS] 
    … …  
    4763                        for f in xrange(len(thickness.weight)): 
    4864 
    49                             self.prg.CoreShellCylinderKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, 
     65                            self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, 
    5066                                    np.float32(axis_theta.value[k]), np.float32(axis_phi.value[l]), np.float32(thickness.value[f]), 
    5167                                    np.float32(length.value[j]), np.float32(radius.value[i]), np.float32(pars['scale']), 
    … …  
    5470                                    np.float32(pars['shell_sld']), np.float32(pars['solvent_sld']),np.uint32(size), 
    5571                                    np.uint32(self.qx.size)) 
    56                             cl.enqueue_copy(self.queue, self.res, self.res_b) 
     72                            cl.enqueue_copy(queue, self.res, self.res_b) 
    5773 
    5874                            sum += self.res 

  • cylcode.py

    r5378e40 r8a20be5  
    66import pyopencl as cl 
    77from weights import GaussianDispersion 
     8from sasmodel import card 
    89 
     10def set_precision(src, qx, qy, dtype): 
     11    qx = np.ascontiguousarray(qx, dtype=dtype) 
     12    qy = np.ascontiguousarray(qy, dtype=dtype) 
     13    if np.dtype(dtype) == np.dtype('float32'): 
     14        header = """\ 
     15#define real float 
     16""" 
     17    else: 
     18        header = """\ 
     19#pragma OPENCL EXTENSION cl_khr_fp64: enable 
     20#define real double 
     21""" 
     22    return header+src, qx, qy 
    923 
    1024class GpuCylinder(object): 
    … …  
    1529    PD_PARS = ['radius', 'length', 'cyl_theta', 'cyl_phi'] 
    1630 
    17     def __init__(self, qx, qy): 
     31    def __init__(self, qx, qy, dtype='float32'): 
    1832 
    19         self.qx = np.asarray(qx, np.float32) 
    20         self.qy = np.asarray(qy, np.float32) 
    2133        #create context, queue, and build program 
    22         self.ctx = cl.create_some_context() 
    23         self.queue = cl.CommandQueue(self.ctx) 
    24         self.prg = cl.Program(self.ctx, open('Kernel-Cylinder.cpp').read()).build() 
     34        ctx,_queue = card() 
     35        src, qx, qy = set_precision(open('NR_BessJ1.cpp').read()+"\n"+open('Kernel-Cylinder.cpp').read(), qx, qy, dtype=dtype) 
     36        self.prg = cl.Program(ctx, src).build() 
     37        self.qx, self.qy = qx, qy 
    2538 
    2639        #buffers 
    2740        mf = cl.mem_flags 
    28         self.qx_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx) 
    29         self.qy_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy) 
    30         self.res_b = cl.Buffer(self.ctx, mf.WRITE_ONLY, qx.nbytes) 
     41        self.qx_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx) 
     42        self.qy_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy) 
     43        self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, qx.nbytes) 
    3144        self.res = np.empty_like(self.qx) 
    3245 
    3346    def eval(self, pars): 
    3447 
    35         radius,length,cyl_theta,cyl_phi = \ 
     48        _ctx,queue = card() 
     49        radius, length, cyl_theta, cyl_phi = \ 
    3650            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) 
    3751             for base in GpuCylinder.PD_PARS] 
    … …  
    4862        sub = pars['sldCyl'] - pars['sldSolv'] 
    4963 
     64        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64 
    5065        #Loop over radius, length, theta, phi weight points 
    5166        for i in xrange(len(radius.weight)): 
    … …  
    5469                    for l in xrange(len(cyl_phi.weight)): 
    5570 
    56                         self.prg.CylinderKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, np.float32(sub), 
    57                                            np.float32(radius.value[i]), np.float32(length.value[j]), np.float32(pars['scale']), 
    58                                            np.float32(radius.weight[i]), np.float32(length.weight[j]), np.float32(cyl_theta.weight[k]), 
    59                                            np.float32(cyl_phi.weight[l]), np.float32(cyl_theta.value[k]), np.float32(cyl_phi.value[l]), 
     71 
     72                        self.prg.CylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, real(sub), 
     73                                           real(radius.value[i]), real(length.value[j]), real(pars['scale']), 
     74                                           real(radius.weight[i]), real(length.weight[j]), real(cyl_theta.weight[k]), 
     75                                           real(cyl_phi.weight[l]), real(cyl_theta.value[k]), real(cyl_phi.value[l]), 
    6076                                           np.uint32(self.qx.size), np.uint32(size)) 
    61                         cl.enqueue_copy(self.queue, self.res, self.res_b) 
     77                        cl.enqueue_copy(queue, self.res, self.res_b) 
    6278                        sum += self.res 
    6379                        vol += radius.weight[i]*length.weight[j]*pow(radius.value[i], 2)*length.value[j] 

  • fit.py

    r5378e40 r8a20be5  
    1313set_beam_stop(data, 0.004) 
    1414 
    15 """ 
     15 
     16 
     17 
    1618model = SasModel(data, GpuCylinder, scale=1, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=0, 
    1719                              cyl_theta=0, cyl_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,length_pd=0.1, 
    1820                              length_pd_n=5, length_pd_nsigma=3, cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
    19                               cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3,) 
    20  
     21                              cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3, dtype='float32') 
     22model.radius.range(0,100) 
     23model.length.range(0, 1000) 
     24model.cyl_theta.range(0,90) 
     25model.cyl_phi.range(0,90) 
     26""" 
    2127model = SasModel(data, GpuEllipse, scale=.027, radius_a=60, radius_b=180, sldEll=.297e-6, sldSolv=5.773e-6, background=4.9, 
    2228                 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, 
    2329                 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, 
    24                  axis_phi_pd_n=6, axis_phi_pd_nsigma=3) 
     30                 axis_phi_pd_n=6, axis_phi_pd_nsigma=3, dtype='float') 
    2531 
    2632model = SasModel(data, GpuLamellar, scale=1, bi_thick=100, sld_bi=.291e-6, sld_sol=5.77e-6, background=0, 
    27                  bi_thick_pd=0.1, bi_thick_pd_n=35, bi_thick_pd_nsigma=3) 
     33                 bi_thick_pd=0.1, bi_thick_pd_n=35, bi_thick_pd_nsigma=3, dtype='float') 
    2834 
    29 """ 
    30 model = SasModel(data, GpuCoreShellCylinder, scale=1, radius=64.1, thickness=1, length=266.96, core_sld=.251e-6, shell_sld=6.2e-6, 
    31                  solvent_sld=5.77e-6, background=0, axis_theta=0, axis_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
     35 
     36model = SasModel(data, GpuCoreShellCylinder, scale=1, radius=64.1, thickness=1, length=266.96, core_sld=1e-6, shell_sld=1e-6, 
     37                 solvent_sld=4e-6, background=0, axis_theta=0, axis_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
    3238                 length_pd=0.1, length_pd_n=10, length_pd_nsigma=3, thickness_pd=0.1, thickness_pd_n=2, thickness_pd_nsigma=3, 
    3339                 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, 
    34                  axis_phi_pd_nsigma=3) 
     40                 axis_phi_pd_nsigma=3, dtype='float') 
     41""" 
    3542 
    3643model.scale.range(0,10) 

  • lamellarcode.py

    r5378e40 r8a20be5  
    33 
    44import numpy as np 
    5 import math 
    65import pyopencl as cl 
    76from weights import GaussianDispersion 
     7 
     8def 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#pragma OPENCL EXTENSION cl_khr_fp64: enable 
     14#define real double 
     15""" 
     16        return header+src,qx,qy 
     17    else: 
     18        return src,qx,qy 
    819 
    920 
    … …  
    1223        'scale':1, 'bi_thick':1, 'sld_bi':1e-6, 'sld_sol':0, 'background':0, 
    1324    } 
    14     PD_PARS = ['bi_thick'] 
    1525 
    16     def __init__(self, qx, qy): 
     26    def __init__(self, qx, qy, dtype='float'): 
    1727 
    18         self.qx = np.asarray(qx, np.float32) 
    19         self.qy = np.asarray(qy, np.float32) 
    2028        #create context, queue, and build program 
    2129        self.ctx = cl.create_some_context() 
    2230        self.queue = cl.CommandQueue(self.ctx) 
    23         self.prg = cl.Program(self.ctx, open('Kernel-Lamellar.cpp').read()).build() 
     31        src,qx,qy = set_precision(open('Kernel-Lamellar.cpp').read(), qx, qy, dtype=dtype) 
     32        self.prg = cl.Program(self.ctx, src).build() 
     33        self.qx, self.qy = qx, qy 
    2434 
    2535        #buffers 
    … …  
    3848        sub = pars['sld_bi'] - pars['sld_sol'] 
    3949 
     50        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64 
    4051        for i in xrange(len(bi_thick.weight)): 
    41             self.prg.LamellarKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, np.float32(bi_thick.value[i]), 
    42                                     np.float32(pars['scale']), np.float32(sub), np.float32(pars['background']), np.uint32(self.qx.size)) 
     52            self.prg.LamellarKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, real(bi_thick.value[i]), 
     53                                    real(pars['scale']), real(sub), real(pars['background']), np.uint32(self.qx.size)) 
    4354            cl.enqueue_copy(self.queue, self.res, self.res_b) 
    4455 
    … …  
    4758 
    4859        return sum/norm + pars['background'] 
    49  
    50     def lamellar_fit(self, pars, b_n=10, b_w=.1, sigma=3): 
    51  
    52         bi_thick = GaussianDispersion(b_n, b_w, sigma) 
    53         bi_thick.value, bi_thick.weight = bi_thick.get_weights(pars.bi_thick, 0, 1000, True) 
    54  
    55         sum, norm = 0.0, 0.0 
    56  
    57         for i in xrange(len(bi_thick.weight)): 
    58             self.prg.LamellarKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, np.float32(bi_thick.value[i]), 
    59                                     np.float32(pars.scale), np.float32(pars.sld_bi), np.float32(pars.sld_sol), 
    60                                     np.float32(pars.background), np.uint32(self.qx.size)) 
    61             cl.enqueue_copy(self.queue, self.res, self.res_b) 
    62  
    63             sum += bi_thick.weight[i]*self.res 
    64             norm += bi_thick.weight[i] 
    65  
    66         return sum/norm + pars.background 
    67  
    6860 
    6961def demo(): 

  • sasmodel.py

    r5378e40 r8a20be5  
    33 
    44import numpy as np 
     5import pyopencl as cl 
    56from bumps.names import Parameter 
    67from sans.dataloader.loader import Loader 
    … …  
    1112    loader = Loader() 
    1213    data = loader.load(filename) 
     14    if data is None: 
     15        raise IOError("Data %r could not be loaded"%filename) 
    1316    return data 
    1417 
    … …  
    4750 
    4851 
     52GPU_CONTEXT = None 
     53GPU_QUEUE = None 
     54def card(): 
     55    global GPU_CONTEXT, GPU_QUEUE 
     56    if GPU_CONTEXT is None: 
     57        GPU_CONTEXT = cl.create_some_context() 
     58        GPU_QUEUE = cl.CommandQueue(GPU_CONTEXT) 
     59    return GPU_CONTEXT, GPU_QUEUE 
     60 
     61 
    4962class SasModel(object): 
    50     def __init__(self, data, model, **kw): 
     63    def __init__(self, data, model, dtype='float32', **kw): 
    5164        self.index = data.mask==0 
    5265        self.iq = data.data[self.index] 
    … …  
    5568        self.qx = data.qx_data 
    5669        self.qy = data.qy_data 
    57         self.gpu = model(self.qx, self.qy) 
     70        self.gpu = model(self.qx, self.qy, dtype=dtype) 
    5871        pd_pars = set(base+attr for base in model.PD_PARS for attr in ('_pd','_pd_n','_pd_nsigma')) 
    5972        total_pars = set(model.PARS.keys()) | pd_pars 
    … …  
    7992    def theory(self): 
    8093        pars = dict((k,v.value) for k,v in self._parameters.items()) 
     94        print pars 
    8195        result = self.gpu.eval(pars) 
    8296        return result 
Note: See TracChangeset for help on using the changeset viewer.