1 | # -*- coding: utf-8 -*- |
---|
2 | ''' |
---|
3 | *GSASIImpsubs - routines used in multiprocessing* |
---|
4 | ------------------------------------------------- |
---|
5 | |
---|
6 | The routines here are called either directly when GSAS-II is used without multiprocessing |
---|
7 | or in separate cores when multiprocessing is used. |
---|
8 | |
---|
9 | These routines are designed to be used in one of two ways: |
---|
10 | |
---|
11 | * when multiprocessing is |
---|
12 | enabled (see global variable useMP) the computational routines are called in |
---|
13 | separate Python interpreter that is created and then deleted after use. |
---|
14 | |
---|
15 | * when useMP is False, these routines are called directly from the main "thread". |
---|
16 | |
---|
17 | Note that :func:`GSASIImpsubs.InitMP` should be called before any of the other routines |
---|
18 | in this module are used. |
---|
19 | ''' |
---|
20 | ########### SVN repository information ################### |
---|
21 | # $Date: $ |
---|
22 | # $Author: $ |
---|
23 | # $Revision: $ |
---|
24 | # $URL: $ |
---|
25 | # $Id: $ |
---|
26 | ########### SVN repository information ################### |
---|
27 | import multiprocessing as mp |
---|
28 | import numpy as np |
---|
29 | import numpy.ma as ma |
---|
30 | import numpy.linalg as nl |
---|
31 | import GSASIIpath |
---|
32 | GSASIIpath.SetVersionNumber("$Revision: 2895 $") |
---|
33 | import GSASIIpwd as G2pwd |
---|
34 | import GSASIIstrMath as G2stMth |
---|
35 | |
---|
36 | sind = lambda x: np.sin(x*np.pi/180.) |
---|
37 | cosd = lambda x: np.cos(x*np.pi/180.) |
---|
38 | tand = lambda x: np.tan(x*np.pi/180.) |
---|
39 | #asind = lambda x: 180.*np.arcsin(x)/np.pi |
---|
40 | #acosd = lambda x: 180.*np.arccos(x)/np.pi |
---|
41 | #atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi |
---|
42 | |
---|
43 | ncores = None |
---|
44 | |
---|
45 | def InitMP(allowMP=True): |
---|
46 | '''Called in routines to initialize use of Multiprocessing |
---|
47 | ''' |
---|
48 | global useMP,ncores |
---|
49 | if ncores is not None: return |
---|
50 | useMP = False |
---|
51 | if not allowMP: |
---|
52 | print('Multiprocessing disabled') |
---|
53 | ncores = 0 |
---|
54 | return |
---|
55 | ncores = GSASIIpath.GetConfigValue('Multiprocessing_cores',-1) |
---|
56 | if ncores < 0: ncores = mp.cpu_count() |
---|
57 | if ncores > 1: |
---|
58 | useMP = True |
---|
59 | #if GSASIIpath.GetConfigValue('debug'): |
---|
60 | if True: |
---|
61 | print('Multiprocessing with {} cores enabled'.format(ncores)) |
---|
62 | |
---|
63 | ################################################################################ |
---|
64 | # derivative computation |
---|
65 | ################################################################################ |
---|
66 | def InitDerivGlobals(im1,calcControls1,SGMT1,hfx1,phfx1,pfx1,G1,GB1,g1,SGData1, |
---|
67 | parmDict1,wave1,shl1,x1,cw1,Ka21,A1,varylist1,dependentVars1, |
---|
68 | dFdvDict1,lamRatio1,kRatio1,doPawley1,pawleyLookup1): |
---|
69 | '''Initialize for the computation of derivatives. Puts lots of junk into the global |
---|
70 | namespace in this module, including the arrays for derivatives (when needed.) |
---|
71 | ''' |
---|
72 | global im,calcControls,SGMT,hfx,phfx,pfx,G,GB,g,SGData,parmDict,wave,shl,x,cw,Ka2,A |
---|
73 | global varylist,dependentVars,dFdvDict,lamRatio,kRatio,doPawley,pawleyLookup |
---|
74 | im = im1 |
---|
75 | calcControls = calcControls1 |
---|
76 | SGMT = SGMT1 |
---|
77 | hfx = hfx1 |
---|
78 | phfx = phfx1 |
---|
79 | pfx = pfx1 |
---|
80 | G = G1 |
---|
81 | GB = GB1 |
---|
82 | g = g1 |
---|
83 | SGData = SGData1 |
---|
84 | parmDict = parmDict1 |
---|
85 | wave = wave1 |
---|
86 | shl = shl1 |
---|
87 | x = ma.getdata(x1) |
---|
88 | cw = cw1 |
---|
89 | Ka2 = Ka21 |
---|
90 | A = A1 |
---|
91 | varylist = varylist1 |
---|
92 | dependentVars = dependentVars1 |
---|
93 | dFdvDict = dFdvDict1 |
---|
94 | lamRatio = lamRatio1 |
---|
95 | kRatio = kRatio1 |
---|
96 | doPawley = doPawley1 |
---|
97 | pawleyLookup = pawleyLookup1 |
---|
98 | # determine the parameters that will have derivatives computed only at end |
---|
99 | global nonatomvarylist |
---|
100 | nonatomvarylist = [] |
---|
101 | for name in varylist: |
---|
102 | if '::RBV;' not in name: |
---|
103 | try: |
---|
104 | aname = name.split(pfx)[1][:2] |
---|
105 | if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc', \ |
---|
106 | 'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param |
---|
107 | except IndexError: |
---|
108 | continue |
---|
109 | nonatomvarylist.append(name) |
---|
110 | global nonatomdependentVars |
---|
111 | nonatomdependentVars = [] |
---|
112 | for name in dependentVars: |
---|
113 | if '::RBV;' not in name: |
---|
114 | try: |
---|
115 | aname = name.split(pfx)[1][:2] |
---|
116 | if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc', \ |
---|
117 | 'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param |
---|
118 | except IndexError: |
---|
119 | continue |
---|
120 | nonatomdependentVars.append(name) |
---|
121 | # create local copies of derivative arrays, if multiprocessing will be used |
---|
122 | if useMP: |
---|
123 | global dMdv |
---|
124 | dMdv = np.zeros(shape=(len(varylist),len(x))) |
---|
125 | global depDerivDict |
---|
126 | depDerivDict = {j:np.zeros(shape=len(x)) for j in dependentVars} |
---|
127 | |
---|
128 | |
---|
129 | def cellVaryDerv(pfx,SGData,dpdA): |
---|
130 | if SGData['SGLaue'] in ['-1',]: |
---|
131 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]], |
---|
132 | [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]] |
---|
133 | elif SGData['SGLaue'] in ['2/m',]: |
---|
134 | if SGData['SGUniq'] == 'a': |
---|
135 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]] |
---|
136 | elif SGData['SGUniq'] == 'b': |
---|
137 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]] |
---|
138 | else: |
---|
139 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]] |
---|
140 | elif SGData['SGLaue'] in ['mmm',]: |
---|
141 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]] |
---|
142 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
143 | return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]] |
---|
144 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
145 | return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]] |
---|
146 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
147 | return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]] |
---|
148 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
149 | return [[pfx+'A0',dpdA[0]]] |
---|
150 | |
---|
151 | def ComputeDerivMPbatch(reflsList): |
---|
152 | '''Computes a the derivatives for a batch of reflections and sums the them into |
---|
153 | global arrays dMdv & depDerivDict. These arrays are returned once the computation |
---|
154 | is completed. |
---|
155 | ''' |
---|
156 | for refl,iref,fmin,fmax,iBeg,iFin in reflsList: |
---|
157 | if ComputeDeriv(refl,iref,fmin,fmax,iBeg,iFin,dMdv,depDerivDict): break |
---|
158 | return dMdv,depDerivDict |
---|
159 | |
---|
160 | def ComputeDeriv(refl,iref,fmin,fmax,iBeg,iFin,dMdv,depDerivDict): |
---|
161 | '''Compute the parameter derivatives for a single reflection and add the results |
---|
162 | into either array dMdv or depDerivDict |
---|
163 | ''' |
---|
164 | global wave |
---|
165 | if im: |
---|
166 | h,k,l,m = refl[:4] |
---|
167 | else: |
---|
168 | h,k,l = refl[:3] |
---|
169 | Uniq = np.inner(refl[:3],SGMT) |
---|
170 | if 'T' in calcControls[hfx+'histType']: |
---|
171 | wave = refl[14+im] |
---|
172 | dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = G2stMth.GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) |
---|
173 | pos = refl[5+im] |
---|
174 | calcKa2 = False |
---|
175 | if 'C' in calcControls[hfx+'histType']: |
---|
176 | tanth = tand(pos/2.0) |
---|
177 | costh = cosd(pos/2.0) |
---|
178 | if Ka2: |
---|
179 | pos2 = refl[5+im]+lamRatio*tanth # + 360/pi * Dlam/lam * tan(th) |
---|
180 | iBeg2 = np.searchsorted(x,pos2-fmin) |
---|
181 | iFin2 = np.searchsorted(x,pos2+fmax) |
---|
182 | if iBeg2-iFin2: |
---|
183 | calcKa2 = True |
---|
184 | iFin = iFin2 |
---|
185 | lenBF = iFin-iBeg |
---|
186 | dMdpk = np.zeros(shape=(6,lenBF)) |
---|
187 | dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin]) |
---|
188 | for i in range(5): |
---|
189 | dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i] |
---|
190 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])} |
---|
191 | if calcKa2: |
---|
192 | dMdpk2 = np.zeros(shape=(6,lenBF)) |
---|
193 | dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,x[iBeg:iFin]) |
---|
194 | for i in range(5): |
---|
195 | dMdpk2[i] = 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i] |
---|
196 | dMdpk2[5] = 100.*cw[iBeg:iFin]*refl[11+im]*dMdipk2[0] |
---|
197 | dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]} |
---|
198 | else: #'T'OF |
---|
199 | lenBF = iFin-iBeg |
---|
200 | if lenBF < 0: return True #bad peak coeff |
---|
201 | dMdpk = np.zeros(shape=(6,lenBF)) |
---|
202 | dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin]) |
---|
203 | for i in range(6): |
---|
204 | dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i] #cw[iBeg:iFin]* |
---|
205 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]} |
---|
206 | if doPawley: |
---|
207 | try: |
---|
208 | if im: |
---|
209 | pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)]) |
---|
210 | else: |
---|
211 | pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]) |
---|
212 | idx = varylist.index(pIdx) |
---|
213 | dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9+im] |
---|
214 | if Ka2: #not for TOF either |
---|
215 | dMdv[idx][iBeg:iFin] += dervDict2['int']/refl[9+im] |
---|
216 | except: # ValueError: |
---|
217 | pass |
---|
218 | if 'C' in calcControls[hfx+'histType']: |
---|
219 | dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = G2stMth.GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict) |
---|
220 | names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'], |
---|
221 | hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'], |
---|
222 | hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'], |
---|
223 | hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'], |
---|
224 | hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'], |
---|
225 | hfx+'DisplaceY':[dpdY,'pos'],} |
---|
226 | if 'Bragg' in calcControls[hfx+'instType']: |
---|
227 | names.update({hfx+'SurfRoughA':[dFdAb[0],'int'], |
---|
228 | hfx+'SurfRoughB':[dFdAb[1],'int'],}) |
---|
229 | else: |
---|
230 | names.update({hfx+'Absorption':[dFdAb,'int'],}) |
---|
231 | else: #'T'OF |
---|
232 | dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = G2stMth.GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict) |
---|
233 | names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'], |
---|
234 | hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'], |
---|
235 | hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'], |
---|
236 | hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'], |
---|
237 | hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'], |
---|
238 | hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'], |
---|
239 | hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],} |
---|
240 | for name in names: |
---|
241 | item = names[name] |
---|
242 | if name in varylist: |
---|
243 | dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]] |
---|
244 | if calcKa2: |
---|
245 | dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict2[item[1]] |
---|
246 | elif name in dependentVars: |
---|
247 | depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]] |
---|
248 | if calcKa2: |
---|
249 | depDerivDict[name][iBeg:iFin] += item[0]*dervDict2[item[1]] |
---|
250 | for iPO in dIdPO: |
---|
251 | if iPO in varylist: |
---|
252 | dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] |
---|
253 | if calcKa2: |
---|
254 | dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict2['int'] |
---|
255 | elif iPO in dependentVars: |
---|
256 | depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] |
---|
257 | if calcKa2: |
---|
258 | depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict2['int'] |
---|
259 | for i,name in enumerate(['omega','chi','phi']): |
---|
260 | aname = pfx+'SH '+name |
---|
261 | if aname in varylist: |
---|
262 | dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int'] |
---|
263 | if calcKa2: |
---|
264 | dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict2['int'] |
---|
265 | elif aname in dependentVars: |
---|
266 | depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int'] |
---|
267 | if calcKa2: |
---|
268 | depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict2['int'] |
---|
269 | for iSH in dFdODF: |
---|
270 | if iSH in varylist: |
---|
271 | dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] |
---|
272 | if calcKa2: |
---|
273 | dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict2['int'] |
---|
274 | elif iSH in dependentVars: |
---|
275 | depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] |
---|
276 | if calcKa2: |
---|
277 | depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict2['int'] |
---|
278 | cellDervNames = cellVaryDerv(pfx,SGData,dpdA) |
---|
279 | for name,dpdA in cellDervNames: |
---|
280 | if name in varylist: |
---|
281 | dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos'] |
---|
282 | if calcKa2: |
---|
283 | dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict2['pos'] |
---|
284 | elif name in dependentVars: #need to scale for mixed phase constraints? |
---|
285 | depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos'] |
---|
286 | if calcKa2: |
---|
287 | depDerivDict[name][iBeg:iFin] += dpdA*dervDict2['pos'] |
---|
288 | dDijDict = G2stMth.GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict) |
---|
289 | for name in dDijDict: |
---|
290 | if name in varylist: |
---|
291 | dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos'] |
---|
292 | if calcKa2: |
---|
293 | dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict2['pos'] |
---|
294 | elif name in dependentVars: |
---|
295 | depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos'] |
---|
296 | if calcKa2: |
---|
297 | depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict2['pos'] |
---|
298 | for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']): |
---|
299 | if name in varylist: |
---|
300 | dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos'] |
---|
301 | if calcKa2: |
---|
302 | dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict2['pos'] |
---|
303 | elif name in dependentVars: |
---|
304 | depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos'] |
---|
305 | if calcKa2: |
---|
306 | depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict2['pos'] |
---|
307 | if 'C' in calcControls[hfx+'histType']: |
---|
308 | sigDict,gamDict = G2stMth.GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict) |
---|
309 | else: #'T'OF |
---|
310 | sigDict,gamDict = G2stMth.GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict) |
---|
311 | for name in gamDict: |
---|
312 | if name in varylist: |
---|
313 | dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam'] |
---|
314 | if calcKa2: |
---|
315 | dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict2['gam'] |
---|
316 | elif name in dependentVars: |
---|
317 | depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam'] |
---|
318 | if calcKa2: |
---|
319 | depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict2['gam'] |
---|
320 | for name in sigDict: |
---|
321 | if name in varylist: |
---|
322 | dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig'] |
---|
323 | if calcKa2: |
---|
324 | dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict2['sig'] |
---|
325 | elif name in dependentVars: |
---|
326 | depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig'] |
---|
327 | if calcKa2: |
---|
328 | depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict2['sig'] |
---|
329 | for name in ['BabA','BabU']: |
---|
330 | if refl[9+im]: |
---|
331 | if phfx+name in varylist: |
---|
332 | dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im] |
---|
333 | if calcKa2: |
---|
334 | dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im] |
---|
335 | elif phfx+name in dependentVars: |
---|
336 | depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im] |
---|
337 | if calcKa2: |
---|
338 | depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im] |
---|
339 | if not doPawley and not parmDict[phfx+'LeBail']: |
---|
340 | #do atom derivatives - for RB,F,X & U so far - how do I scale mixed phase constraints? |
---|
341 | corr = 0. |
---|
342 | #corr2 = 0. |
---|
343 | if refl[9+im]: |
---|
344 | corr = dervDict['int']/refl[9+im] |
---|
345 | #if calcKa2: # commented out in Bob's code. Why? |
---|
346 | # corr2 = dervDict2['int']/refl[9+im] |
---|
347 | for name in nonatomvarylist: |
---|
348 | dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr |
---|
349 | #if calcKa2: |
---|
350 | # dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr2 # unneeded w/o above |
---|
351 | for name in nonatomdependentVars: |
---|
352 | depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr |
---|
353 | #if calcKa2: |
---|
354 | # depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr2 |
---|
355 | |
---|
356 | ################################################################################ |
---|
357 | # Fobs Squared computation |
---|
358 | ################################################################################ |
---|
359 | #x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2 |
---|
360 | def InitFobsSqGlobals(x1,ratio1,shl1,xB1,xF1,im1,lamRatio1,kRatio1,xMask1,Ka21): |
---|
361 | '''Initialize for the computation of Fobs Squared for powder histograms. |
---|
362 | Puts lots of junk into the global namespace in this module. |
---|
363 | ''' |
---|
364 | global x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2 |
---|
365 | x = ma.getdata(x1) |
---|
366 | ratio = ratio1 |
---|
367 | shl = shl1 |
---|
368 | xB = xB1 |
---|
369 | xF = xF1 |
---|
370 | im = im1 |
---|
371 | lamRatio = lamRatio1 |
---|
372 | kRatio = kRatio1 |
---|
373 | xMask = xMask1 |
---|
374 | Ka2 = Ka21 |
---|
375 | |
---|
376 | def ComputeFobsSqCWbatch(profList): |
---|
377 | sInt = 0 |
---|
378 | resList = [] |
---|
379 | for refl,iref in profList: |
---|
380 | icod = ComputeFobsSqCW(refl,iref) |
---|
381 | if type(icod) is tuple: |
---|
382 | resList.append((icod[0],iref)) |
---|
383 | sInt += icod[1] |
---|
384 | elif icod == -1: |
---|
385 | res.append((None,iref)) |
---|
386 | elif icod == -2: |
---|
387 | break |
---|
388 | return sInt,resList |
---|
389 | |
---|
390 | def ComputeFobsSqTOFbatch(profList): |
---|
391 | sInt = 0 |
---|
392 | resList = [] |
---|
393 | for refl,iref in profList: |
---|
394 | icod = ComputeFobsSqTOF(refl,iref) |
---|
395 | if type(icod) is tuple: |
---|
396 | resList.append((icod[0],iref)) |
---|
397 | sInt += icod[1] |
---|
398 | elif icod == -1: |
---|
399 | res.append((None,iref)) |
---|
400 | elif icod == -2: |
---|
401 | break |
---|
402 | return sInt,resList |
---|
403 | |
---|
404 | def ComputeFobsSqCW(refl,iref): |
---|
405 | yp = np.zeros(len(x)) # not masked |
---|
406 | sInt = 0 |
---|
407 | refl8im = 0 |
---|
408 | Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl) |
---|
409 | iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin)) |
---|
410 | iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF)) |
---|
411 | iFin2 = iFin |
---|
412 | if not iBeg+iFin: #peak below low limit - skip peak |
---|
413 | return 0 |
---|
414 | if ma.all(xMask[iBeg:iFin]): #peak entirely masked - skip peak |
---|
415 | return -1 |
---|
416 | elif not iBeg-iFin: #peak above high limit - done |
---|
417 | return -2 |
---|
418 | elif iBeg < iFin: |
---|
419 | yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3( |
---|
420 | refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin]) |
---|
421 | sInt = refl[11+im]*refl[9+im] |
---|
422 | if Ka2: |
---|
423 | pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0) # + 360/pi * Dlam/lam * tan(th) |
---|
424 | Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl) |
---|
425 | iBeg2 = max(xB,np.searchsorted(x,pos2-fmin)) |
---|
426 | iFin2 = min(np.searchsorted(x,pos2+fmax),xF) |
---|
427 | if iFin2 > iBeg2: |
---|
428 | yp[iBeg2:iFin2] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3( |
---|
429 | pos2,refl[6+im],refl[7+im],shl,x[iBeg2:iFin2]) |
---|
430 | sInt *= 1.+kRatio |
---|
431 | refl8im = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0)) |
---|
432 | return refl8im,sInt |
---|
433 | |
---|
434 | def ComputeFobsSqTOF(refl,iref): |
---|
435 | yp = np.zeros(len(x)) # not masked |
---|
436 | refl8im = 0 |
---|
437 | Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im]) |
---|
438 | iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin)) |
---|
439 | iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF)) |
---|
440 | if not iBeg+iFin: #peak below low limit - skip peak |
---|
441 | return 0 |
---|
442 | if ma.all(xMask[iBeg:iFin]): #peak entirely masked - skip peak |
---|
443 | return -1 |
---|
444 | elif not iBeg-iFin: #peak above high limit - done |
---|
445 | return -2 |
---|
446 | if iBeg < iFin: |
---|
447 | yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt( |
---|
448 | refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin]) |
---|
449 | refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0)) |
---|
450 | return refl8im,refl[11+im]*refl[9+im] |
---|
451 | ################################################################################ |
---|
452 | # Powder Profile computation |
---|
453 | ################################################################################ |
---|
454 | def InitPwdrProfGlobals(im1,shl1,x1): |
---|
455 | '''Initialize for the computation of Fobs Squared for powder histograms. |
---|
456 | Puts lots of junk into the global namespace in this module. |
---|
457 | ''' |
---|
458 | global im,shl,x |
---|
459 | im = im1 |
---|
460 | shl = shl1 |
---|
461 | x = ma.getdata(x1) |
---|
462 | global cw |
---|
463 | cw = np.diff(x) |
---|
464 | cw = np.append(cw,cw[-1]) |
---|
465 | # create local copies of ycalc array |
---|
466 | if useMP: |
---|
467 | global yc |
---|
468 | yc = np.zeros_like(x1) |
---|
469 | |
---|
470 | |
---|
471 | def ComputePwdrProfCW(profList): |
---|
472 | 'Compute the peaks profile for a set of CW peaks and add into the yc array' |
---|
473 | for pos,refl,iBeg,iFin,kRatio in profList: |
---|
474 | yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3( |
---|
475 | pos,refl[6+im],refl[7+im],shl,x[iBeg:iFin]) |
---|
476 | return yc |
---|
477 | |
---|
478 | def ComputePwdrProfTOF(profList): |
---|
479 | 'Compute the peaks profile for a set of TOF peaks and add into the yc array' |
---|
480 | for pos,refl,iBeg,iFin in profList: |
---|
481 | yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt( |
---|
482 | pos,refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])/cw[iBeg:iFin] |
---|
483 | return yc |
---|
484 | |
---|