source: sasmodels/sasmodel.py @ ce27e21

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since ce27e21 was ce27e21, checked in by Paul Kienzle <pkienzle@…>, 10 years ago

first pass for sasview wrapper around opencl models

  • Property mode set to 100644
File size: 7.2 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3
4import datetime
5import numpy as np
6
7
8TIC = None
9def tic():
10    global TIC
11    TIC = datetime.datetime.now()
12
13def toc():
14    now = datetime.datetime.now()
15    return (now-TIC).total_seconds()
16
17def load_data(filename):
18    from sans.dataloader.loader import Loader
19    loader = Loader()
20    data = loader.load(filename)
21    if data is None:
22        raise IOError("Data %r could not be loaded"%filename)
23    return data
24
25def set_precision(src, qx, qy, dtype):
26    qx = np.ascontiguousarray(qx, dtype=dtype)
27    qy = np.ascontiguousarray(qy, dtype=dtype)
28    if np.dtype(dtype) == np.dtype('float32'):
29        header = """\
30
31#define REAL(x) (x##f)
32#define real float
33"""
34    else:
35        header = """\
36#pragma OPENCL EXTENSION cl_khr_fp64: enable
37#define REAL(x) (x)
38#define real double
39"""
40    return header+src, qx, qy
41
42def set_precision_1d(src, q, dtype):
43    q = np.ascontiguousarray(q, dtype=dtype)
44    if np.dtype(dtype) == np.dtype('float32'):
45        header = """\
46#define real float
47"""
48    else:
49        header = """\
50#pragma OPENCL EXTENSION cl_khr_fp64: enable
51#define real double
52"""
53    return header+src, q
54
55def set_beam_stop(data, radius, outer=None):
56    from sans.dataloader.manipulations import Ringcut
57    if hasattr(data, 'qx_data'):
58        data.mask = Ringcut(0, radius)(data)
59        if outer is not None:
60            data.mask += Ringcut(outer,np.inf)(data)
61    else:
62        data.mask = (data.x>=radius)
63        if outer is not None:
64            data.mask &= (data.x<outer)
65
66def set_half(data, half):
67    from sans.dataloader.manipulations import Boxcut
68    if half == 'right':
69        data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data)
70    if half == 'left':
71        data.mask += Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data)
72
73def set_top(data, max):
74    from sans.dataloader.manipulations import Boxcut
75    data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data)
76
77def plot_data(data, iq, vmin=None, vmax=None):
78    from numpy.ma import masked_array
79    import matplotlib.pyplot as plt
80    img = masked_array(iq, data.mask)
81    xmin, xmax = min(data.qx_data), max(data.qx_data)
82    ymin, ymax = min(data.qy_data), max(data.qy_data)
83    plt.imshow(img.reshape(128,128),
84               interpolation='nearest', aspect=1, origin='upper',
85               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
86
87def plot_result2D(data, theory, view='linear'):
88    import matplotlib.pyplot as plt
89    from numpy.ma import masked_array, masked
90    #print "not a number",sum(np.isnan(data.data))
91    #data.data[data.data<0.05] = 0.5
92    mdata = masked_array(data.data, data.mask)
93    mdata[np.isnan(mdata)] = masked
94    if view is 'log':
95        mdata[mdata <= 0] = masked
96        mdata = np.log10(mdata)
97        mtheory = masked_array(np.log10(theory), mdata.mask)
98    else:
99        mtheory = masked_array(theory, mdata.mask)
100    mresid = masked_array((theory-data.data)/data.err_data, data.mask)
101    vmin = min(mdata.min(), mtheory.min())
102    vmax = max(mdata.max(), mtheory.max())
103    print np.exp(np.mean(mtheory)), np.std(mtheory),np.max(mtheory),np.min(mtheory)
104
105    plt.subplot(1, 3, 1)
106    plot_data(data, mdata, vmin=vmin, vmax=vmax)
107    plt.colorbar()
108    plt.subplot(1, 3, 2)
109    plot_data(data, mtheory, vmin=vmin, vmax=vmax)
110    plt.colorbar()
111    plt.subplot(1, 3, 3)
112    print abs(mresid).max()
113    plot_data(data, mresid)
114    plt.colorbar()
115
116
117def plot_result1D(data, theory, view='linear'):
118    import matplotlib.pyplot as plt
119    from numpy.ma import masked_array, masked
120    #print "not a number",sum(np.isnan(data.y))
121    #data.y[data.y<0.05] = 0.5
122    mdata = masked_array(data.y, data.mask)
123    mdata[np.isnan(mdata)] = masked
124    if view is 'log':
125        mdata[mdata <= 0] = masked
126    mtheory = masked_array(theory, mdata.mask)
127    mresid = masked_array((theory-data.y)/data.dy, mdata.mask)
128
129    plt.subplot(1,2,1)
130    plt.errorbar(data.x, mdata, yerr=data.dy)
131    plt.plot(data.x, mtheory, '-', hold=True)
132    plt.yscale(view)
133    plt.subplot(1, 2, 2)
134    plt.plot(data.x, mresid, 'x')
135    #plt.axhline(1, color='black', ls='--',lw=1, hold=True)
136    #plt.axhline(0, color='black', lw=1, hold=True)
137    #plt.axhline(-1, color='black', ls='--',lw=1, hold=True)
138
139
140
141GPU_CONTEXT = None
142GPU_QUEUE = None
143def card():
144    global GPU_CONTEXT, GPU_QUEUE
145    if GPU_CONTEXT is None:
146        import pyopencl as cl
147        GPU_CONTEXT = cl.create_some_context()
148        GPU_QUEUE = cl.CommandQueue(GPU_CONTEXT)
149    return GPU_CONTEXT, GPU_QUEUE
150
151
152class SasModel(object):
153    def __init__(self, data, model, dtype='float32', **kw):
154        from bumps.names import Parameter
155
156        self.__dict__['_parameters'] = {}
157        #self.name = data.filename
158        self.is2D = hasattr(data,'qx_data')
159        self.data = data
160        if self.is2D:
161            self.index = (data.mask==0) & (~np.isnan(data.data))
162            self.iq = data.data[self.index]
163            self.diq = data.err_data[self.index]
164            self.qx = data.qx_data
165            self.qy = data.qy_data
166            self.gpu = model(self.qx, self.qy, dtype=dtype)
167        else:
168            self.index = (data.mask==0) & (~np.isnan(data.y))
169            self.iq = data.y[self.index]
170            self.diq = data.dy[self.index]
171            self.q = data.x
172            self.gpu = model(self.q, dtype=dtype)
173        pd_pars = set(base+attr for base in model.PD_PARS for attr in ('_pd','_pd_n','_pd_nsigma'))
174        total_pars = set(model.PARS.keys()) | pd_pars
175        extra_pars = set(kw.keys()) - total_pars
176        if extra_pars:
177            raise TypeError("unexpected parameters %s"%(str(extra_pars,)))
178        pars = model.PARS.copy()
179        pars.update((base+'_pd', 0) for base in model.PD_PARS)
180        pars.update((base+'_pd_n', 35) for base in model.PD_PARS)
181        pars.update((base+'_pd_nsigma', 3) for base in model.PD_PARS)
182        pars.update(kw)
183        for k,v in pars.items():
184            setattr(self, k, Parameter.default(v, name=k))
185        self._parameter_names = set(pars.keys())
186        self.update()
187
188    def update(self):
189        self._cache = {}
190
191    def numpoints(self):
192        return len(self.iq)
193
194    def parameters(self):
195        return dict((k,getattr(self,k)) for k in self._parameter_names)
196
197    def theory(self):
198        if 'theory' not in self._cache:
199            pars = dict((k,getattr(self,k).value) for k in self._parameter_names)
200            #print pars
201            self._cache['theory'] = self.gpu.eval(pars)
202        return self._cache['theory']
203
204    def residuals(self):
205        #if np.any(self.err ==0): print "zeros in err"
206        return (self.theory()[self.index]-self.iq)/self.diq
207
208    def nllf(self):
209        R = self.residuals()
210        #if np.any(np.isnan(R)): print "NaN in residuals"
211        return 0.5*np.sum(R**2)
212
213    def __call__(self):
214        return 2*self.nllf()/self.dof
215
216    def plot(self, view='log'):
217        if self.is2D:
218            plot_result2D(self.data, self.theory(), view=view)
219        else:
220            plot_result1D(self.data, self.theory(), view=view)
221
222    def save(self, basename):
223        pass
224
225def demo():
226    data = load_data('DEC07086.DAT')
227    set_beam_stop(data, 0.004)
228    plot_data(data)
229    import matplotlib.pyplot as plt; plt.show()
230
231if __name__ == "__main__":
232    demo()
Note: See TracBrowser for help on using the repository browser.