source: branch/2frame/GSASIImpsubs.py @ 2941

Last change on this file since 2941 was 2941, checked in by toby, 6 years ago

prepare for multiprocessing implementation
  • Property svn:eol-style set to native
File size: 16.2 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIImpsubs - routines used in multiprocessing*
4-------------------------------------------------
5
6The routines here are called either directly when GSAS-II is used without multiprocessing
7or in separate cores when multiprocessing is used.
8
9These 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
17Note that :func:`GSASIImpsubs.InitMP` should be called before any of the other routines
18in this module are used.
19'''
20########### SVN repository information ###################
21# $Date: $
22# $Author: $
23# $Revision: $
24# $URL: $
25# $Id: $
26########### SVN repository information ###################
27import multiprocessing as mp
28import numpy as np
29import numpy.ma as ma
30import GSASIIpath
31GSASIIpath.SetVersionNumber("$Revision: 2895 $")
32import GSASIIpwd as G2pwd
33import GSASIIstrMath as G2stMth
34
35cosd = lambda x: np.cos(x*np.pi/180.)
36tand = lambda x: np.tan(x*np.pi/180.)
37#acosd = lambda x: 180.*np.arccos(x)/np.pi
38#atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
39   
40ncores = None
41
42def InitMP():
43    '''Called in routines to initialize use of Multiprocessing
44    '''
45    global useMP,ncores
46    if ncores is not None: return
47    useMP = False
48    ncores = GSASIIpath.GetConfigValue('Multiprocessing_cores',-1)
49    if ncores < 0: ncores = mp.cpu_count()
50    if ncores > 1:
51        useMP = True
52    #if GSASIIpath.GetConfigValue('debug'):
53    if True:
54        print('Multiprocessing with {} cores enabled'.format(ncores))
55
56################################################################################
57# derivative computation
58################################################################################       
59def InitDerivGlobals(im1,calcControls1,SGMT1,hfx1,phfx1,pfx1,G1,GB1,g1,SGData1,
60                     parmDict1,wave1,shl1,x1,cw1,Ka21,A1,varylist1,dependentVars1,
61                     dFdvDict1,lamRatio1,kRatio1,doPawley1):
62    '''Initialize for the computation of derivatives. Puts lots of junk into the global
63    namespace in this module, including the arrays for derivatives (when needed.)
64    '''
65    global im,calcControls,SGMT,hfx,phfx,pfx,G,GB,g,SGData,parmDict,wave,shl,x,cw,Ka2,A
66    global varylist,dependentVars,dFdvDict,lamRatio,kRatio,doPawley
67    im = im1
68    calcControls = calcControls1
69    SGMT = SGMT1
70    hfx = hfx1
71    phfx = phfx1
72    pfx = pfx1
73    G = G1
74    GB = GB1
75    g = g1
76    SGData = SGData1
77    parmDict = parmDict1
78    wave = wave1
79    shl = shl1
80    x = ma.getdata(x1)
81    cw = cw1
82    Ka2 = Ka21
83    A = A1
84    varylist = varylist1
85    dependentVars = dependentVars1
86    dFdvDict = dFdvDict1
87    lamRatio = lamRatio1
88    kRatio = kRatio1
89    doPawley = doPawley1
90    # determine the parameters that will have derivatives computed only at end
91    global nonatomvarylist
92    nonatomvarylist = []
93    for name in varylist:
94        if '::RBV;' not in name:
95            try:
96                aname = name.split(pfx)[1][:2]
97                if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
98                    'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
99            except IndexError:
100                continue
101        nonatomvarylist.append(name)
102    global nonatomdependentVars
103    nonatomdependentVars = []
104    for name in dependentVars:
105        if '::RBV;' not in name:
106            try:
107                aname = name.split(pfx)[1][:2]
108                if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
109                    'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
110            except IndexError:
111                continue
112        nonatomdependentVars.append(name)
113    # create local copies of derivative arrays, if multiprocessing will be used
114    if useMP:
115        global dMdv
116        dMdv = np.zeros(shape=(len(varylist),len(x)))
117        global depDerivDict
118        depDerivDict = {j:np.zeros(shape=len(x)) for j in dependentVars}
119           
120
121def cellVaryDerv(pfx,SGData,dpdA): 
122    if SGData['SGLaue'] in ['-1',]:
123        return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
124            [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
125    elif SGData['SGLaue'] in ['2/m',]:
126        if SGData['SGUniq'] == 'a':
127            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
128        elif SGData['SGUniq'] == 'b':
129            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
130        else:
131            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
132    elif SGData['SGLaue'] in ['mmm',]:
133        return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
134    elif SGData['SGLaue'] in ['4/m','4/mmm']:
135        return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
136    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
137        return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
138    elif SGData['SGLaue'] in ['3R', '3mR']:
139        return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
140    elif SGData['SGLaue'] in ['m3m','m3']:
141        return [[pfx+'A0',dpdA[0]]]
142
143def ComputeDerivMPbatch(reflsList):
144    '''Computes a the derivatives for a batch of reflections and sums the them into
145    global arrays dMdv & depDerivDict. These arrays are returned once the computation
146    is completed.
147    '''
148    for refl,iref,fmin,fmax,iBeg,iFin in reflsList:
149        if ComputeDeriv(refl,iref,fmin,fmax,iBeg,iFin,dMdv,depDerivDict): break
150    return dMdv,depDerivDict
151
152def ComputeDeriv(refl,iref,fmin,fmax,iBeg,iFin,dMdv,depDerivDict):
153    '''Compute the parameter derivatives for a single reflection and add the results
154    into either array dMdv or depDerivDict
155    '''
156    global wave
157    if im:
158        h,k,l,m = refl[:4]
159    else:
160        h,k,l = refl[:3]
161    Uniq = np.inner(refl[:3],SGMT)
162    if 'T' in calcControls[hfx+'histType']:
163        wave = refl[14+im]
164    dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = G2stMth.GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
165    pos = refl[5+im]
166    if 'C' in calcControls[hfx+'histType']:
167        tanth = tand(pos/2.0)
168        costh = cosd(pos/2.0)
169        calcKa2 = False
170        if Ka2:
171            pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
172            iBeg2 = np.searchsorted(x,pos2-fmin)
173            iFin2 = np.searchsorted(x,pos2+fmax)
174            if iBeg2-iFin2:
175                calcKa2 = True
176                iFin = iFin2       
177        lenBF = iFin-iBeg
178        dMdpk = np.zeros(shape=(6,lenBF))
179        dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin])
180        for i in range(5):
181            dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
182        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
183        if calcKa2: 
184            dMdpk2 = np.zeros(shape=(6,lenBF))
185            dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,x[iBeg:iFin])
186            for i in range(5):
187                dMdpk2[i] = 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
188            dMdpk2[5] = 100.*cw[iBeg:iFin]*refl[11+im]*dMdipk2[0]
189            dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
190    else:   #'T'OF
191        lenBF = iFin-iBeg
192        if lenBF < 0: return True  #bad peak coeff
193        dMdpk = np.zeros(shape=(6,lenBF))
194        dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])
195        for i in range(6):
196            dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
197        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
198    if doPawley:
199        try:
200            if im:
201                pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
202            else:
203                pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
204            idx = varylist.index(pIdx)
205            dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9+im]
206            if Ka2: #not for TOF either
207                dMdv[idx][iBeg:iFin] += dervDict2['int']/refl[9+im]
208        except: # ValueError:
209            pass
210    if 'C' in calcControls[hfx+'histType']:
211        dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = G2stMth.GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
212        names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
213            hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
214            hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
215            hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
216            hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
217            hfx+'DisplaceY':[dpdY,'pos'],}
218        if 'Bragg' in calcControls[hfx+'instType']:
219            names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
220                hfx+'SurfRoughB':[dFdAb[1],'int'],})
221        else:
222            names.update({hfx+'Absorption':[dFdAb,'int'],})
223    else:   #'T'OF
224        dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
225        names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
226            hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
227            hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
228            hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
229            hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
230            hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
231            hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
232    for name in names:
233        item = names[name]
234        if name in varylist:
235            dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
236            if calcKa2:
237                dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict2[item[1]]
238        elif name in dependentVars:
239            depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
240            if calcKa2:
241                depDerivDict[name][iBeg:iFin] += item[0]*dervDict2[item[1]]
242    for iPO in dIdPO:
243        if iPO in varylist:
244            dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
245            if calcKa2:
246                dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict2['int']
247        elif iPO in dependentVars:
248            depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
249            if calcKa2:
250                depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict2['int']
251    for i,name in enumerate(['omega','chi','phi']):
252        aname = pfx+'SH '+name
253        if aname in varylist:
254            dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
255            if calcKa2:
256                dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict2['int']
257        elif aname in dependentVars:
258            depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
259            if calcKa2:
260                depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict2['int']
261    for iSH in dFdODF:
262        if iSH in varylist:
263            dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
264            if calcKa2:
265                dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict2['int']
266        elif iSH in dependentVars:
267            depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
268            if calcKa2:
269                depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict2['int']
270    cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
271    for name,dpdA in cellDervNames:
272        if name in varylist:
273            dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
274            if calcKa2:
275                dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict2['pos']
276        elif name in dependentVars: #need to scale for mixed phase constraints?
277            depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
278            if calcKa2:
279                depDerivDict[name][iBeg:iFin] += dpdA*dervDict2['pos']
280    dDijDict = G2stMth.GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
281    for name in dDijDict:
282        if name in varylist:
283            dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
284            if calcKa2:
285                dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict2['pos']
286        elif name in dependentVars:
287            depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
288            if calcKa2:
289                depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict2['pos']
290    for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
291        if name in varylist:
292            dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
293            if calcKa2:
294                dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict2['pos']
295        elif name in dependentVars:
296            depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
297            if calcKa2:
298                depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict2['pos']
299    if 'C' in calcControls[hfx+'histType']:
300        sigDict,gamDict = G2stMth.GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
301    else:   #'T'OF
302        sigDict,gamDict = G2stMth.GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
303    for name in gamDict:
304        if name in varylist:
305            dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
306            if calcKa2:
307                dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict2['gam']
308        elif name in dependentVars:
309            depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
310            if calcKa2:
311                depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict2['gam']
312    for name in sigDict:
313        if name in varylist:
314            dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
315            if calcKa2:
316                dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict2['sig']
317        elif name in dependentVars:
318            depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
319            if calcKa2:
320                depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict2['sig']
321    for name in ['BabA','BabU']:
322        if refl[9+im]:
323            if phfx+name in varylist:
324                dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
325                if calcKa2:
326                    dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]
327            elif phfx+name in dependentVars:                   
328                depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
329                if calcKa2:
330                    depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]                 
331    if not doPawley and not parmDict[phfx+'LeBail']:
332        #do atom derivatives -  for RB,F,X & U so far - how do I scale mixed phase constraints?
333        corr = 0.
334        #corr2 = 0.
335        if refl[9+im]:             
336            corr = dervDict['int']/refl[9+im]
337            #if calcKa2:  # commented out in Bob's code. Why?
338            #    corr2 = dervDict2['int']/refl[9+im]
339        for name in nonatomvarylist:
340            dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
341            #if calcKa2:
342            #   dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr2 # unneeded w/o above
343        for name in nonatomdependentVars:
344           depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
345           #if calcKa2:
346           #    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr2
Note: See TracBrowser for help on using the repository browser.