Changeset d772f5d in sasmodels


Ignore:
Timestamp:
Jul 23, 2014 1:25:44 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:
dbb0048
Parents:
9a9f5b5
Message:

Added 1D Fit, fixed fitting error

Files:
3 added
1 deleted
11 edited

Legend:

Unmodified
Added
Removed

  • Capcyl_Kfun.cpp

    r2de9a5e rd772f5d  
    1414        return(val); 
    1515} 
    16 /* 
    17 real Gauss76Z[76]={ 
    18          .999505948362153*(-1.0),               //0 
    19          .997397786355355*(-1.0), 
    20          .993608772723527*(-1.0), 
    21          .988144453359837*(-1.0), 
    22          .981013938975656*(-1.0), 
    23          .972229228520377*(-1.0), 
    24          .961805126758768*(-1.0), 
    25          .949759207710896*(-1.0), 
    26          .936111781934811*(-1.0), 
    27          .92088586125215*(-1.0), 
    28          .904107119545567*(-1.0),               //10 
    29          .885803849292083*(-1.0), 
    30          .866006913771982*(-1.0), 
    31          .844749694983342*(-1.0), 
    32          .822068037328975*(-1.0), 
    33          .7980001871612*(-1.0), 
    34          .77258672828181*(-1.0), 
    35          .74587051350361*(-1.0), 
    36          .717896592387704*(-1.0), 
    37          .688712135277641*(-1.0), 
    38          .658366353758143*(-1.0),               //20 
    39          .626910417672267*(-1.0), 
    40          .594397368836793*(-1.0), 
    41          .560882031601237*(-1.0), 
    42          .526420920401243*(-1.0), 
    43          .491072144462194*(-1.0), 
    44          .454895309813726*(-1.0), 
    45          .417951418780327*(-1.0), 
    46          .380302767117504*(-1.0), 
    47          .342012838966962*(-1.0), 
    48          .303146199807908*(-1.0),               //30 
    49          .263768387584994*(-1.0), 
    50          .223945802196474*(-1.0), 
    51          .183745593528914*(-1.0), 
    52          .143235548227268*(-1.0), 
    53          .102483975391227*(-1.0), 
    54          .0615595913906112*(-1.0), 
    55          .0205314039939986*(-1.0), 
    56         -.0205314039939986*(-1.0), 
    57         -.0615595913906112*(-1.0), 
    58         -.102483975391227*(-1.0),                       //40 
    59         -.143235548227268*(-1.0), 
    60         -.183745593528914*(-1.0), 
    61         -.223945802196474*(-1.0), 
    62         -.263768387584994*(-1.0), 
    63         -.303146199807908*(-1.0), 
    64         -.342012838966962*(-1.0), 
    65         -.380302767117504*(-1.0), 
    66         -.417951418780327*(-1.0), 
    67         -.454895309813726*(-1.0), 
    68         -.491072144462194*(-1.0),               //50 
    69         -.526420920401243*(-1.0), 
    70         -.560882031601237*(-1.0), 
    71         -.594397368836793*(-1.0), 
    72         -.626910417672267*(-1.0), 
    73         -.658366353758143*(-1.0), 
    74         -.688712135277641*(-1.0), 
    75         -.717896592387704*(-1.0), 
    76         -.74587051350361*(-1.0), 
    77         -.77258672828181*(-1.0), 
    78         -.7980001871612*(-1.0), //60 
    79         -.822068037328975*(-1.0), 
    80         -.844749694983342*(-1.0), 
    81         -.866006913771982*(-1.0), 
    82         -.885803849292083*(-1.0), 
    83         -.904107119545567*(-1.0), 
    84         -.92088586125215*(-1.0), 
    85         -.936111781934811*(-1.0), 
    86         -.949759207710896*(-1.0), 
    87         -.961805126758768*(-1.0), 
    88         -.972229228520377*(-1.0),               //70 
    89         -.981013938975656*(-1.0), 
    90         -.988144453359837*(-1.0), 
    91         -.993608772723527*(-1.0), 
    92         -.997397786355355*(-1.0), 
    93         -.999505948362153*(-1.0)                //75 
    94 }; 
    95 real Gauss76Wt[76]={ 
    96         .00126779163408536,             //0 
    97         .00294910295364247, 
    98         .00462793522803742, 
    99         .00629918049732845, 
    100         .00795984747723973, 
    101         .00960710541471375, 
    102         .0112381685696677, 
    103         .0128502838475101, 
    104         .0144407317482767, 
    105         .0160068299122486, 
    106         .0175459372914742,              //10 
    107         .0190554584671906, 
    108         .020532847967908, 
    109         .0219756145344162, 
    110         .0233813253070112, 
    111         .0247476099206597, 
    112         .026072164497986, 
    113         .0273527555318275, 
    114         .028587223650054, 
    115         .029773487255905, 
    116         .0309095460374916,              //20 
    117         .0319934843404216, 
    118         .0330234743977917, 
    119         .0339977794120564, 
    120         .0349147564835508, 
    121         .0357728593807139, 
    122         .0365706411473296, 
    123         .0373067565423816, 
    124         .0379799643084053, 
    125         .0385891292645067, 
    126         .0391332242205184,              //30 
    127         .0396113317090621, 
    128         .0400226455325968, 
    129         .040366472122844, 
    130         .0406422317102947, 
    131         .0408494593018285, 
    132         .040987805464794, 
    133         .0410570369162294, 
    134         .0410570369162294, 
    135         .040987805464794, 
    136         .0408494593018285,              //40 
    137         .0406422317102947, 
    138         .040366472122844, 
    139         .0400226455325968, 
    140         .0396113317090621, 
    141         .0391332242205184, 
    142         .0385891292645067, 
    143         .0379799643084053, 
    144         .0373067565423816, 
    145         .0365706411473296, 
    146         .0357728593807139,              //50 
    147         .0349147564835508, 
    148         .0339977794120564, 
    149         .0330234743977917, 
    150         .0319934843404216, 
    151         .0309095460374916, 
    152         .029773487255905, 
    153         .028587223650054, 
    154         .0273527555318275, 
    155         .026072164497986, 
    156         .0247476099206597,              //60 
    157         .0233813253070112, 
    158         .0219756145344162, 
    159         .020532847967908, 
    160         .0190554584671906, 
    161         .0175459372914742, 
    162         .0160068299122486, 
    163         .0144407317482767, 
    164         .0128502838475101, 
    165         .0112381685696677, 
    166         .00960710541471375,             //70 
    167         .00795984747723973, 
    168         .00629918049732845, 
    169         .00462793522803742, 
    170         .00294910295364247, 
    171         .00126779163408536              //75 (indexed from 0) 
    172 }; 
    173 */ 

  • Kernel-CapCyl.cpp

    r09e15be rd772f5d  
    1616        real cos_val = cyl_x*qx[i]/q + cyl_y*qy[i]/q; 
    1717        real alpha = acos(cos_val); 
    18         real yyy=0; real ans1=0; real ans2=0; real y=0; real xx=0; real ans=0; real zij=0; real be=0; real summj=0; 
     18        real ans1=0; real ans2=0; real y=0; real xx=0; real ans=0; real zij=0; real be=0; real summj=0; 
    1919 
    2020        real hDist = -1.0*sqrt(fabs(rad_cap*rad_cap-rad_cyl*rad_cyl)); 
    … …  
    2525        { 
    2626            zij = (Gauss76Z[j]*(1.0-vaj)+vaj+1.0)/2.0; 
    27             yyy = Gauss76Wt[j]*ConvLens_kernel(length,rad_cyl,rad_cap,q,zij,alpha); 
    28             summj += yyy; 
     27            summj += Gauss76Wt[j]*ConvLens_kernel(length,rad_cyl,rad_cap,q,zij,alpha); 
    2928        } 
    30         real inner = (1.0-vaj)/2.0*summj*4.0*pi*rad_cap*rad_cap*rad_cap; 
     29        real yyy = (1.0-vaj)/2.0*summj*4.0*pi*rad_cap*rad_cap*rad_cap; 
    3130        real arg1 = q*length/2.0*cos(alpha); 
    3231        real arg2 = q*rad_cyl*sin(alpha); 
    33         yyy = inner; 
    3432 
    3533        if(arg2 == 0) {be = 0.5;} 

  • code_capcyl.py

    r8a21ba3 rd772f5d  
    77from weights import GaussianDispersion 
    88from sasmodel import card 
    9 from Capcyl_Gauss import Gauss76Wt, Gauss76Z 
     9from Gauss import Gauss76Wt, Gauss76Z 
    1010 
    1111def set_precision(src, qx, qy, dtype): 
    … …  
    4343        self.qx_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx) 
    4444        self.qy_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy) 
    45         G, Z = Gauss76Wt, Gauss76Z 
     45        G = np.ascontiguousarray(Gauss76Wt, dtype=dtype) 
     46        Z = np.ascontiguousarray(Gauss76Z, dtype=dtype) 
    4647        self.Gauss76W_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=G) 
    4748        self.Gauss76Z_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=Z) 

  • code_coreshellcyl.py

    r09e15be rd772f5d  
    4848                                     for base in GpuCoreShellCylinder.PD_PARS] 
    4949 
    50         radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 1000, True) 
    51         length.value, length.weight = length.get_weights(pars['length'], 0, 1000, True) 
    52         thickness.value, thickness.weight = thickness.get_weights(pars['thickness'], 0, 1000, True) 
     50        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True) 
     51        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True) 
     52        thickness.value, thickness.weight = thickness.get_weights(pars['thickness'], 0, 10000, True) 
    5353        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False) 
    5454        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False) 
    5555 
    5656        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0 
     57 
     58        print radius.value 
     59        print thickness.weight 
     60        print axis_phi.weight 
     61        print axis_theta.weight 
     62        print length.value 
     63 
    5764        size = len(axis_theta.weight) 
    5865 
    … …  
    8087                                    *axis_phi.weight[p] 
    8188 
    82         if size>1: 
    83             norm /= math.asin(1.0) 
     89        #if size>1: 
     90         #   norm /= math.asin(1.0) 
    8491        if vol != 0.0 and norm_vol != 0.0: 
    8592            sum *= norm_vol/vol 

  • code_cylinder.py

    r09e15be rd772f5d  
    33 
    44import numpy as np 
    5 import math 
    65import pyopencl as cl 
    76from weights import GaussianDispersion 
    … …  
    2120""" 
    2221    return header+src, qx, qy 
     22 
     23def set_precision_1d(src, q, dtype): 
     24    q = np.ascontiguousarray(q, dtype=dtype) 
     25    if np.dtype(dtype) == np.dtype('float32'): 
     26        header = """\ 
     27#define real float 
     28""" 
     29    else: 
     30        header = """\ 
     31#pragma OPENCL EXTENSION cl_khr_fp64: enable 
     32#define real double 
     33""" 
     34    return header+src, q 
    2335 
    2436class GpuCylinder(object): 
    … …  
    5264 
    5365        #Get the weights for each 
    54         radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 1000, True) 
    55         length.value, length.weight = length.get_weights(pars['length'], 0, 1000, True) 
     66        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True) 
     67        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True) 
    5668        cyl_theta.value, cyl_theta.weight = cyl_theta.get_weights(pars['cyl_theta'], -90, 180, False) 
    5769        cyl_phi.value, cyl_phi.weight = cyl_phi.get_weights(pars['cyl_phi'], -90, 180, False) 
    … …  
    8698        return sum/norm+pars['background'] 
    8799 
    88 def demo(): 
    89     from time import time 
    90     import matplotlib.pyplot as plt 
     100class OneDGpuCylinder(object): 
     101    PARS = { 
     102        'scale':1,'radius':1,'length':1,'sldCyl':1e-6,'sldSolv':0,'background':0, 
     103        'bolim':0, 'uplim':90 
     104    } 
     105    PD_PARS = ['radius', 'length'] 
    91106 
    92     #create qx and qy evenly spaces 
    93     qx = np.linspace(-.02, .02, 128) 
    94     qy = np.linspace(-.02, .02, 128) 
    95     qx, qy = np.meshgrid(qx, qy) 
     107    def __init__(self, q, dtype='float32'): 
    96108 
    97     #saved shape of qx 
    98     r_shape = qx.shape 
     109        #create context, queue, and build program 
     110        ctx,_queue = card() 
     111        trala = open('NR_BessJ1.cpp').read()+"\n"+open('OneDCyl_Kfun.cpp').read()+"\n"+open('Kernel-OneDCylinder.cpp').read() 
     112        src, self.q = set_precision_1d(trala, q, dtype=dtype) 
     113        self.prg = cl.Program(ctx, src).build() 
    99114 
    100     #reshape for calculation; resize as float32 
    101     qx = qx.flatten() 
    102     qy = qy.flatten() 
     115        #buffers 
     116        mf = cl.mem_flags 
     117        self.q_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.q) 
     118        self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, q.nbytes) 
     119        self.res = np.empty_like(self.q) 
    103120 
    104     pars = CylinderParameters(scale=1, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=0, 
    105                               cyl_theta=0, cyl_phi=0) 
    106     t = time() 
    107     result = GpuCylinder(qx, qy) 
    108     result.x = result.cylinder_fit(pars, r_n=10, t_n=10, l_n=10, p_n=10, r_w=.1, t_w=.1, l_w=.1, p_w=.1, sigma=3) 
    109     result.x = np.reshape(result.x, r_shape) 
    110     tt = time() 
     121    def eval(self, pars): 
    111122 
    112     print("Time taken: %f" % (tt - t)) 
     123        _ctx,queue = card() 
     124        radius, length = \ 
     125            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) 
     126             for base in OneDGpuCylinder.PD_PARS] 
    113127 
    114     plt.pcolormesh(result.x) 
    115     plt.show() 
     128        #Get the weights for each 
     129        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True) 
     130        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True) 
    116131 
     132        #Perform the computation, with all weight points 
     133        sum, norm, vol = 0.0, 0.0, 0.0, 
     134        sub = pars['sldCyl'] - pars['sldSolv'] 
    117135 
    118 if __name__=="__main__": 
    119     demo() 
     136        real = np.float32 if self.q.dtype == np.dtype('float32') else np.float64 
     137        #Loop over radius, length, theta, phi weight points 
     138        for r in xrange(len(radius.weight)): 
     139            for l in xrange(len(length.weight)): 
     140                        self.prg.OneDCylKernel(queue, self.q.shape, None, self.q_b, self.res_b, real(sub), 
     141                                           real(length.value[l]), real(radius.value[r]), real(pars['scale']), 
     142                                           np.uint32(self.q.size), real(pars['uplim']), real(pars['bolim'])) 
     143                        cl.enqueue_copy(queue, self.res, self.res_b) 
     144                        sum += radius.weight[r]*length.weight[l]*self.res*pow(radius.value[r],2)*length.value[l] 
     145                        vol += radius.weight[r]*length.weight[l] *pow(radius.value[r],2)*length.value[l] 
     146                        norm += radius.weight[r]*length.weight[l] 
     147 
     148        if vol != 0.0 and norm != 0.0: 
     149            sum *= norm/vol 
     150 
     151        return sum/norm + pars['background'] 
     152 
     153        
     154 
     155         
    120156 
    121157 
    … …  
    145181 
    146182 
     183 

  • code_ellipse.py

    r09e15be rd772f5d  
    4949             for base in GpuEllipse.PD_PARS] 
    5050 
    51         radius_a.value, radius_a.weight = radius_a.get_weights(pars['radius_a'], 0, 1000, True) 
    52         radius_b.value, radius_b.weight = radius_b.get_weights(pars['radius_b'], 0, 1000, True) 
     51        radius_a.value, radius_a.weight = radius_a.get_weights(pars['radius_a'], 0, 10000, True) 
     52        radius_b.value, radius_b.weight = radius_b.get_weights(pars['radius_b'], 0, 10000, True) 
    5353        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False) 
    5454        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False) 

  • code_lamellar.py

    r2de9a5e rd772f5d  
    4545 
    4646        bi_thick = GaussianDispersion(int(pars['bi_thick_pd_n']), pars['bi_thick_pd'], pars['bi_thick_pd_nsigma']) 
    47         bi_thick.value, bi_thick.weight = bi_thick.get_weights(pars['bi_thick'], 0, 1000, True) 
     47        bi_thick.value, bi_thick.weight = bi_thick.get_weights(pars['bi_thick'], 0, 10000, True) 
    4848 
    4949        sum, norm = 0.0, 0.0 

  • code_triaxialellipse.py

    r09e15be rd772f5d  
    4848             for base in GpuTriEllipse.PD_PARS] 
    4949 
    50         axisA.value, axisA.weight = axisA.get_weights(pars['axisA'], 0, 1000, True) 
    51         axisB.value, axisB.weight = axisB.get_weights(pars['axisB'], 0, 1000, True) 
    52         axisC.value, axisC.weight = axisC.get_weights(pars['axisC'], 0, 1000, True) 
     50        axisA.value, axisA.weight = axisA.get_weights(pars['axisA'], 0, 10000, True) 
     51        axisB.value, axisB.weight = axisB.get_weights(pars['axisB'], 0, 10000, True) 
     52        axisC.value, axisC.weight = axisC.get_weights(pars['axisC'], 0, 10000, True) 
    5353        theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False) 
    5454        phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False) 

  • fit.py

    r79fcc40 rd772f5d  
    33 
    44from bumps.names import * 
    5 from code_cylinder import GpuCylinder 
     5from code_cylinder import GpuCylinder, OneDGpuCylinder 
    66from code_lamellar import GpuLamellar 
    77from code_ellipse import GpuEllipse 
    … …  
    1010from code_triaxialellipse import GpuTriEllipse 
    1111from sasmodel import SasModel, load_data, set_beam_stop, set_half 
     12import numpy as np 
     13 
     14""" IMPORT THE DATA USED """ 
     15 
     16data = load_data('December/Tangential/Sector0/DEC07133.ABS') 
     17#data = load_data('December/DEC07133.DAT') 
     18 
     19""" SET INNER BEAM STOP, OUTER RING, AND MASK HALF OF THE DATA """ 
     20set_beam_stop(data, 0.0052)#, outer=0.025) 
     21#set_half(data, 'left') 
    1222 
    1323 
    14 data = load_data('DEC07282.DAT') 
    15 set_beam_stop(data, 0.0048)#, outer=0.025) 
    16 #set_half(data, 'right') 
    17 """ 
    18 model = SasModel(data, GpuCylinder, scale=0.09402, radius=106, length=442.298, 
    19                                     sldCyl=.291e-6, sldSolv=5.77e-6, background=0, 
    20                                     cyl_theta=0, cyl_phi=0, 
    21                                     radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=0, 
    22                                     length_pd=0.1, length_pd_n=5, length_pd_nsigma=0, 
    23                                     cyl_theta_pd=0.1, cyl_theta_pd_n=10, cyl_theta_pd_nsigma=3, 
    24                                     cyl_phi_pd=0.1, cyl_phi_pd_n=4, cyl_phi_pd_nsigma=3, 
    25                                     dtype='float') 
    26 model.radius.range(0, 1000) 
    27 model.length.range(0, 1000) 
    28 #model.cyl_theta.range(0,90) 
    29 #model.cyl_phi.range(0,90) 
    30 model.scale.range(0, 1) 
    31 """ 
     24 
     25model = SasModel(data, OneDGpuCylinder, 
     26scale=0.0013, 
     27radius=105, 
     28length=1000, 
     29background=21, 
     30sldCyl=.291e-6,sldSolv=5.77e-6, 
     31radius_pd=0.1,radius_pd_n=10,radius_pd_nsigma=0, 
     32length_pd=0.1,length_pd_n=5,length_pd_nsigma=0, 
     33bolim=0.0, 
     34uplim=90) #bottom limit, upper limit of angle integral 
     35 
     36 
    3237""" 
    3338model = SasModel(data, GpuEllipse, 
    34                  scale=0.0007, radius_a=42.8251, radius_b=700, sldEll=.291e-6, sldSolv=7.105e-6, 
    35                  background=10, axis_theta=0, axis_phi=0, 
    36                  radius_a_pd=0.222296, radius_a_pd_n=1, radius_a_pd_nsigma=0, 
    37                  radius_b_pd=.000128, radius_b_pd_n=1, radius_b_pd_nsigma=0, 
    38                  axis_theta_pd=24.8059, axis_theta_pd_n=40, axis_theta_pd_nsigma=3, 
    39                  axis_phi_pd=2.63698e-05, axis_phi_pd_n=20, axis_phi_pd_nsigma=3, 
    40                  dtype='float') 
     39scale=0.0011, 
     40radius_a=100, radius_b=800.8, 
     41sldEll=.291e-6, sldSolv=7.105e-6, 
     42background=8.30161, 
     43axis_theta=0, axis_phi=0, 
     44axis_theta_pd=20, axis_theta_pd_n=40, axis_theta_pd_nsigma=3, 
     45radius_a_pd=0.222296, radius_a_pd_n=1, radius_a_pd_nsigma=0, 
     46radius_b_pd=.000128, radius_b_pd_n=1, radius_b_pd_nsigma=0, 
     47axis_phi_pd=2.63698e-05, axis_phi_pd_n=20, axis_phi_pd_nsigma=0, 
     48dtype='float') 
    4149 
     50 
     51# SET THE FITTING PARAMETERS 
    4252model.radius_a.range(15, 1000) 
    43 #model.radius_b.range(15, 1000) 
     53model.radius_b.range(15, 1000) 
    4454#model.axis_theta_pd.range(0, 360) 
    4555#model.background.range(0,1000) 
    4656model.scale.range(0, 1) 
    4757""" 
     58 
    4859""" 
    49 model = SasModel(data, GpuLamellar, scale=0.0128, bi_thick=35.1, sld_bi=.291e-6, sld_sol=5.77e-6, background=17.3, 
    50                  bi_thick_pd= 0.000009, bi_thick_pd_n=35, bi_thick_pd_nsigma=3, dtype='float') 
     60model = SasModel(data, GpuLamellar, 
     61scale=0.70, 
     62bi_thick=5, 
     63sld_bi=.291e-6,sld_sol=5.77e-6, 
     64background=85.23, 
     65bi_thick_pd= 0.0013, bi_thick_pd_n=5, bi_thick_pd_nsigma=3, 
     66dtype='float') 
     67 
     68# SET THE FITTING PARAMETERS 
    5169model.bi_thick.range(0, 1000) 
    5270model.scale.range(0, 1) 
    … …  
    5472#model.background.range(0, 1000) 
    5573""" 
     74 
     75 
     76 
     77""" 
     78model = SasModel(data, GpuCylinder, 
     79scale=0.0013, radius=105, length=1000, 
     80sldCyl=.291e-6, sldSolv=5.77e-6, background=21, 
     81cyl_theta=90, cyl_phi=0, 
     82cyl_theta_pd=534, cyl_theta_pd_n=40, cyl_theta_pd_nsigma=3, 
     83 
     84# SET THE FITTING PARAMETERS 
     85radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=0, 
     86length_pd=0.1, length_pd_n=5, length_pd_nsigma=0, 
     87cyl_phi_pd=0.1, cyl_phi_pd_n=4, cyl_phi_pd_nsigma=0, 
     88dtype='float') 
     89#model.radius.range(0, 1000) 
     90#model.length.range(0, 1000) 
     91#model.cyl_theta_pd.range(0,90) 
     92model.scale.range(0, 1) 
     93model.background.range(0, 1000) 
     94""" 
     95 
    5696""" 
    5797model = SasModel(data, GpuCoreShellCylinder, 
    58                  scale= 1.77881e-06, radius=325, thickness=25, length=34.2709, 
    59                  core_sld=1e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 
    60                  background=223.827, axis_theta=90, axis_phi=0, 
    61                  axis_theta_pd=15.8, 
     98                 scale= 0.08, radius=200, thickness=30, length=2000, 
     99                 core_sld=7e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 
     100                 background=0, axis_theta=0, axis_phi=0, 
    62101 
    63                  radius_pd=0.1, radius_pd_n=1, radius_pd_nsigma=0, 
    64                  length_pd=0.1, length_pd_n=1, length_pd_nsigma=0, 
     102                 radius_pd=0.38, radius_pd_n=10, radius_pd_nsigma=3, 
     103                 length_pd=.9, length_pd_n=10, length_pd_nsigma=3, 
    65104                 thickness_pd=0.1, thickness_pd_n=1, thickness_pd_nsigma=0, 
    66                  axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
    67                  axis_phi_pd=0.0008748, axis_phi_pd_n=10, axis_phi_pd_nsigma=3, 
     105                 axis_theta_pd=10, axis_theta_pd_n=40, axis_theta_pd_nsigma=3, 
     106                 axis_phi_pd=0.1, axis_phi_pd_n=1, axis_phi_pd_nsigma=0, 
    68107                 dtype='float') 
    69108 
     109# SET THE FITTING PARAMETERS 
    70110model.radius.range(15, 1000) 
    71111#model.length.range(0, 1000) 
    … …  
    80120""" 
    81121 
     122""" 
    82123 
    83 """ 
    84124model = SasModel(data, GpuCapCylinder, scale=1, rad_cyl=20, rad_cap=40, length=400, sld_capcyl=1e-6, sld_solv=6.3e-6, 
    85125                 background=0, theta=0, phi=0, rad_cyl_pd=.1, rad_cyl_pd_n=10, rad_cyl_nsigma=3, rad_cap_pd=.1, rad_cap_pd_n=1, 
    … …  
    87127                 theta_pd_nsigma=3, phi_pd=.1, phi_pd_n=4, phi_pd_nsigma=3, dtype='float') 
    88128""" 
    89  
     129""" 
    90130 
    91131model = SasModel(data, GpuTriEllipse, 
    … …  
    100140                 axisC_pd=.1, axisC_pd_n=1, axisC_pd_nsigma=0, dtype='float') 
    101141 
    102 #model.axisA.range(15, 1000) 
    103 #model.axisB.range(15, 1000) 
     142# SET THE FITTING PARAMETERS 
     143model.axisA.range(15, 1000) 
     144model.axisB.range(15, 1000) 
    104145#model.axisC.range(15, 1000) 
    105146#model.background.range(0,1000) 
    … …  
    110151 
    111152 
    112  
    113  
     153""" 
    114154 
    115155problem = FitProblem(model) 
    116156 
    117  

  • multi-fit.py

    r9a9f5b5 rd772f5d  
    55from code_lamellar import GpuLamellar 
    66from code_ellipse import GpuEllipse 
    7 from sasmodel import SasModel, load_data, set_beam_stop, set_half 
     7from multisasmodels import SasModel, load_data, set_beam_stop, set_half 
    88import numpy as np 
    99 
    10 data = load_data('DEC07269.DAT') 
     10data = load_data('DEC07282.DAT') 
    1111set_beam_stop(data, 0.0052)#, outer=0.025) 
    1212#set_half(data, 'left') 
    1313 
    1414truth = SasModel(data, GpuEllipse, 
    15                 scale=.5, radius_a=45.265, radius_b=600.8, sldEll=.291e-6, sldSolv=7.105e-6, 
     15                scale=.0011, radius_a=45.265, radius_b=600.8, sldEll=.291e-6, sldSolv=7.105e-6, 
    1616                 background=8.30161, axis_theta=0, axis_phi=0, 
    1717                 radius_a_pd=0.222296, radius_a_pd_n=1, radius_a_pd_nsigma=0, 

  • sasmodel.py

    r9a9f5b5 rd772f5d  
    1919 
    2020def set_beam_stop(data, radius, outer=None): 
    21     data.mask = Ringcut(0, radius)(data) 
    22     if outer is not None: 
    23         data.mask += Ringcut(outer,np.inf)(data) 
     21    if hasattr(data, 'qx_data'): 
     22        data.mask = Ringcut(0, radius)(data) 
     23        if outer is not None: 
     24            data.mask += Ringcut(outer,np.inf)(data) 
     25    else: 
     26        data.mask = (data.x>=radius) 
     27        if outer is not None: 
     28            data.mask &= (data.x<outer) 
    2429 
    2530def set_half(data, half): 
    … …  
    4146               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
    4247 
    43 def plot_result(data, theory, view='linear'): 
     48def plot_result2D(data, theory, view='linear'): 
    4449    import matplotlib.pyplot as plt 
    4550    from numpy.ma import masked_array, masked 
    46     plt.subplot(1, 3, 1) 
    4751    #print "not a number",sum(np.isnan(data.data)) 
    4852    #data.data[data.data<0.05] = 0.5 
    … …  
    5963    vmax = max(mdata.max(), mtheory.max()) 
    6064 
     65    plt.subplot(1, 3, 1) 
    6166    plot_data(data, mdata, vmin=vmin, vmax=vmax) 
    6267    plt.colorbar() 
    … …  
    6974 
    7075 
    71 def demo(): 
    72     data = load_data('JUN03289.DAT') 
    73     set_beam_stop(data, 0.004) 
    74     plot_data(data) 
    75     import matplotlib.pyplot as plt; plt.show() 
     76def plot_result1D(data, theory, view='linear'): 
     77    import matplotlib.pyplot as plt 
     78    from numpy.ma import masked_array, masked 
     79    #print "not a number",sum(np.isnan(data.y)) 
     80    #data.y[data.y<0.05] = 0.5 
     81    mdata = masked_array(data.y, data.mask) 
     82    mdata[np.isnan(mdata)] = masked 
     83    if view is 'log': 
     84        mdata[mdata <= 0] = masked 
     85    mtheory = masked_array(theory, mdata.mask) 
     86    mresid = masked_array((theory-data.y)/data.dy, mdata.mask) 
     87 
     88    plt.subplot(1,2,1) 
     89    plt.errorbar(data.x, mdata, yerr=data.dy) 
     90    plt.plot(data.x, mtheory, '-', hold=True) 
     91    plt.yscale(view) 
     92    plt.subplot(1, 2, 2) 
     93    plt.plot(data.x, mresid, 'x') 
     94    #plt.axhline(1, color='black', ls='--',lw=1, hold=True) 
     95    #plt.axhline(0, color='black', lw=1, hold=True) 
     96    #plt.axhline(-1, color='black', ls='--',lw=1, hold=True) 
     97 
    7698 
    7799 
    … …  
    89111    def __init__(self, data, model, dtype='float32', **kw): 
    90112        self.__dict__['_parameters'] = {} 
    91         self.index = (data.mask==0) & (~np.isnan(data.data)) 
    92         self.iq = data.data[self.index] 
    93         self.diq = data.err_data[self.index] 
     113        self.is2D = hasattr(data,'qx_data') 
    94114        self.data = data 
    95         self.qx = data.qx_data 
    96         self.qy = data.qy_data 
    97         self.gpu = model(self.qx, self.qy, dtype=dtype) 
     115        if self.is2D: 
     116            self.index = (data.mask==0) & (~np.isnan(data.data)) 
     117            self.iq = data.data[self.index] 
     118            self.diq = data.err_data[self.index] 
     119            self.qx = data.qx_data 
     120            self.qy = data.qy_data 
     121            self.gpu = model(self.qx, self.qy, dtype=dtype) 
     122        else: 
     123            self.index = (data.mask==0) & (~np.isnan(data.y)) 
     124            self.iq = data.y[self.index] 
     125            self.diq = data.dy[self.index] 
     126            self.q = data.x 
     127            self.gpu = model(self.q, dtype=dtype) 
    98128        pd_pars = set(base+attr for base in model.PD_PARS for attr in ('_pd','_pd_n','_pd_nsigma')) 
    99129        total_pars = set(model.PARS.keys()) | pd_pars 
    … …  
    107137        pars.update(kw) 
    108138        self._parameters = dict((k, Parameter.default(v, name=k)) for k, v in pars.items()) 
    109  
    110     def set_result(self, result): 
    111         self.result = result 
    112         return self.result 
    113  
    114     def get_result(self): 
    115         return self.result 
    116139 
    117140    def numpoints(self): 
    … …  
    137160    def residuals(self): 
    138161        #if np.any(self.err ==0): print "zeros in err" 
    139         return (self.get_result()[self.index]-self.iq)/self.diq 
     162        return (self.theory()[self.index]-self.iq)/self.diq 
    140163 
    141164    def nllf(self): 
    … …  
    148171 
    149172    def plot(self, view='log'): 
    150         plot_result(self.data, self.get_result(), view=view) 
     173        if self.is2D: 
     174            plot_result2D(self.data, self.theory(), view=view) 
     175        else: 
     176            plot_result1D(self.data, self.theory(), view=view) 
    151177 
    152178    def save(self, basename): 
    … …  
    155181    def update(self): 
    156182        pass 
     183 
     184 
     185def demo(): 
     186    data = load_data('DEC07086.DAT') 
     187    set_beam_stop(data, 0.004) 
     188    plot_data(data) 
     189    import matplotlib.pyplot as plt; plt.show() 
     190 
     191if __name__ == "__main__": 
     192    demo() 
Note: See TracChangeset for help on using the changeset viewer.