source: trunk/GSASIImpsubs.py @ 5114

Last change on this file since 5114 was 5114, checked in by vondreele, 9 months ago

various fixes for EDX data fitting by LeBail?
various fixes for SAD analysis - plotting updates & importer
removed a couple of unused lambdas from G2sasd

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 10.0 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: 2021-12-21 14:53:00 +0000 (Tue, 21 Dec 2021) $
22# $Author: vondreele $
23# $Revision: 5114 $
24# $URL: trunk/GSASIImpsubs.py $
25# $Id: GSASIImpsubs.py 5114 2021-12-21 14:53:00Z vondreele $
26########### SVN repository information ###################
27from __future__ import division, print_function
28import multiprocessing as mp
29import numpy as np
30import numpy.ma as ma
31import GSASIIpath
32GSASIIpath.SetVersionNumber("$Revision: 5114 $")
33import GSASIIpwd as G2pwd
34import GSASIIfiles as G2fil
35
36sind = lambda x: np.sin(x*np.pi/180.)
37cosd = lambda x: np.cos(x*np.pi/180.)
38tand = 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   
43ncores = None
44
45def ResetMP():
46    '''Call after changing Config var 'Multiprocessing_cores' to force a resetting
47    of the useMP from the parameter.
48    '''
49    global ncores
50    ncores = None
51   
52def InitMP(allowMP=True):
53    '''Called to initialize use of Multiprocessing
54    '''
55    global useMP,ncores
56    if ncores is not None: return useMP,ncores
57    useMP = False
58    if not allowMP:
59        G2fil.G2Print('Multiprocessing disabled')
60        ncores = 0
61        return useMP,ncores
62    ncores = GSASIIpath.GetConfigValue('Multiprocessing_cores',0)
63    if ncores < 0: ncores = mp.cpu_count()//2
64    if ncores > 1:
65        useMP = True
66    if useMP:
67        G2fil.G2Print('Multiprocessing with {} cores enabled'.format(ncores))
68    return useMP,ncores
69
70################################################################################
71# Fobs Squared computation
72################################################################################       
73def InitFobsSqGlobals(x1,ratio1,shl1,xB1,xF1,im1,lamRatio1,kRatio1,xMask1,Ka21):
74    '''Initialize for the computation of Fobs Squared for powder histograms.
75    Puts lots of junk into the global namespace in this module.
76    '''
77    global x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2,cw
78    x = ma.getdata(x1)
79    cw = np.diff(x)
80    cw = np.append(cw,cw[-1])
81    ratio = ratio1
82    shl = shl1
83    xB = xB1
84    xF = xF1
85    im = im1
86    lamRatio = lamRatio1
87    kRatio = kRatio1
88    xMask = xMask1
89    Ka2 = Ka21
90
91def ComputeFobsSqCWbatch(profList):
92    sInt = 0
93    resList = []
94    for refl,iref in profList:
95        icod = ComputeFobsSqCW(refl,iref)
96        if type(icod) is tuple:
97            resList.append((icod[0],iref))
98            sInt += icod[1]
99        elif icod == -1:
100            resList.append((None,iref))
101        elif icod == -2:
102            break
103    return sInt,resList
104
105def ComputeFobsSqTOFbatch(profList):
106    sInt = 0
107    resList = []
108    for refl,iref in profList:
109        icod = ComputeFobsSqTOF(refl,iref)
110        if type(icod) is tuple:
111            resList.append((icod[0],iref))
112            sInt += icod[1]
113        elif icod == -1:
114            resList.append((None,iref))
115        elif icod == -2:
116            break
117    return sInt,resList
118       
119def ComputeFobsSqPinkbatch(profList):
120    sInt = 0
121    resList = []
122    for refl,iref in profList:
123        icod = ComputeFobsSqPink(refl,iref)
124        if type(icod) is tuple:
125            resList.append((icod[0],iref))
126            sInt += icod[1]
127        elif icod == -1:
128            resList.append((None,iref))
129        elif icod == -2:
130            break
131    return sInt,resList
132
133def ComputeFobsSqEDbatch(profList):
134    sInt = 0
135    resList = []
136    for refl,iref in profList:
137        icod = ComputeFobsSqED(refl,iref)
138        if type(icod) is tuple:
139            resList.append((icod[0],iref))
140            sInt += icod[1]
141        elif icod == -1:
142            resList.append((None,iref))
143        elif icod == -2:
144            break
145    return sInt,resList
146
147def ComputeFobsSqCW(refl,iref):
148    yp = np.zeros(len(x)) # not masked
149    sInt = 0
150    refl8im = 0
151    Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
152    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
153    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
154    iFin2 = iFin
155    if not iBeg+iFin:       #peak below low limit - skip peak
156        return 0
157    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
158        return -1
159    elif not iBeg-iFin:     #peak above high limit - done
160        return -2
161    elif iBeg < iFin:
162        fp,sumfp = G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin])
163        yp[iBeg:iFin] = 100.*refl[11+im]*refl[9+im]*fp*cw[iBeg:iFin]/sumfp
164        sInt = refl[11+im]*refl[9+im]
165        if Ka2:
166            pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
167            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
168            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
169            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
170            if iFin2 > iBeg2:
171                fp2,sumfp2 = G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,x[iBeg2:iFin2])
172                yp[iBeg2:iFin2] += 100.*refl[11+im]*refl[9+im]*kRatio*fp2*cw[iBeg2:iFin2]/sumfp2
173                sInt *= 1.+kRatio
174    refl8im = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
175    return refl8im,sInt
176
177def ComputeFobsSqTOF(refl,iref):
178    yp = np.zeros(len(x)) # not masked
179    refl8im = 0
180    Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
181    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
182    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
183    if not iBeg+iFin:       #peak below low limit - skip peak
184        return 0
185    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
186        return -1
187    elif not iBeg-iFin:     #peak above high limit - done
188        return -2
189    if iBeg < iFin:
190        fp,sumfp = G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])
191        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*fp*cw[iBeg:iFin]/sumfp
192    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
193    return refl8im,refl[11+im]*refl[9+im]
194
195def ComputeFobsSqPink(refl,iref):
196    yp = np.zeros(len(x)) # not masked
197    refl8im = 0
198    Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.)
199    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
200    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
201    if not iBeg+iFin:       #peak below low limit - skip peak
202        return 0
203    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
204        return -1
205    elif not iBeg-iFin:     #peak above high limit - done
206        return -2
207    if iBeg < iFin:
208        fp,sumfp = G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.,x[iBeg:iFin])
209        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*fp*cw[iBeg:iFin]/sumfp
210    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
211    return refl8im,refl[11+im]*refl[9+im]
212
213def ComputeFobsSqED(refl,iref):
214    yp = np.zeros(len(x)) # not masked
215    refl8im = 0
216    Wd,fmin,fmax = G2pwd.getWidthsED(refl[5+im],refl[6+im])
217    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
218    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
219    if not iBeg+iFin:       #peak below low limit - skip peak
220        return 0
221    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
222        return -1
223    elif not iBeg-iFin:     #peak above high limit - done
224        return -2
225    if iBeg < iFin:
226        fp,sumfp = G2pwd.getPsVoigt(refl[5+im],refl[6+im]*1.e4,0.001,x[iBeg:iFin])
227        yp[iBeg:iFin] = 100.*refl[9+im]*fp*cw[iBeg:iFin]/sumfp
228    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin],0.0))
229    return refl8im,refl[9+im]
230################################################################################
231# Powder Profile computation
232################################################################################       
233def InitPwdrProfGlobals(im1,shl1,x1):
234    '''Initialize for the computation of Fobs Squared for powder histograms.
235    Puts lots of junk into the global namespace in this module.
236    '''
237    global im,shl,x,cw,yc
238    im = im1
239    shl = shl1
240    x = ma.getdata(x1)
241    cw = np.diff(x)
242    cw = np.append(cw,cw[-1])
243    # create local copies of ycalc array
244    yc = np.zeros_like(x)
245
246def ComputePwdrProfCW(profList):
247    'Compute the peaks profile for a set of CW peaks and add into the yc array'
248    for pos,refl,iBeg,iFin,kRatio in profList:
249        fp = G2pwd.getFCJVoigt3(pos,refl[6+im],refl[7+im],shl,x[iBeg:iFin])[0]
250        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*fp/cw[iBeg:iFin]
251    return yc
252
253def ComputePwdrProfTOF(profList):
254    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
255    for pos,refl,iBeg,iFin in profList:
256        fp = G2pwd.getEpsVoigt(pos,refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])[0]
257        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*fp*cw[iBeg:iFin]
258   
259def ComputePwdrProfPink(profList):
260    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
261    for pos,refl,iBeg,iFin in profList:
262        fp = G2pwd.getEpsVoigt(pos,refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.,x[iBeg:iFin])[0]
263        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*fp
264    return yc
265
266def ComputePwdrProfED(profList):
267    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
268    for pos,refl,iBeg,iFin in profList:
269        fp = G2pwd.getPsVoigt(pos,refl[6+im]*1.e4,0.001,x[iBeg:iFin])[0]
270        yc[iBeg:iFin] += refl[9+im]*fp
271    return yc
Note: See TracBrowser for help on using the repository browser.