source: trunk/GSASIImpsubs.py @ 4538

Last change on this file since 4538 was 4521, checked in by vondreele, 5 years ago

complete pink Rietveld refinement
fix problem with weights
enhance testDeriv
put new pink parms in parm dictionary
fix problem with Instrument parm display after load
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 8.4 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: 2020-07-15 20:08:04 +0000 (Wed, 15 Jul 2020) $
22# $Author: vondreele $
23# $Revision: 4521 $
24# $URL: trunk/GSASIImpsubs.py $
25# $Id: GSASIImpsubs.py 4521 2020-07-15 20:08:04Z 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: 4521 $")
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
78    x = ma.getdata(x1)
79    ratio = ratio1
80    shl = shl1
81    xB = xB1
82    xF = xF1
83    im = im1
84    lamRatio = lamRatio1
85    kRatio = kRatio1
86    xMask = xMask1
87    Ka2 = Ka21
88
89def ComputeFobsSqCWbatch(profList):
90    sInt = 0
91    resList = []
92    for refl,iref in profList:
93        icod = ComputeFobsSqCW(refl,iref)
94        if type(icod) is tuple:
95            resList.append((icod[0],iref))
96            sInt += icod[1]
97        elif icod == -1:
98            resList.append((None,iref))
99        elif icod == -2:
100            break
101    return sInt,resList
102
103def ComputeFobsSqTOFbatch(profList):
104    sInt = 0
105    resList = []
106    for refl,iref in profList:
107        icod = ComputeFobsSqTOF(refl,iref)
108        if type(icod) is tuple:
109            resList.append((icod[0],iref))
110            sInt += icod[1]
111        elif icod == -1:
112            resList.append((None,iref))
113        elif icod == -2:
114            break
115    return sInt,resList
116       
117def ComputeFobsSqPinkbatch(profList):
118    sInt = 0
119    resList = []
120    for refl,iref in profList:
121        icod = ComputeFobsSqPink(refl,iref)
122        if type(icod) is tuple:
123            resList.append((icod[0],iref))
124            sInt += icod[1]
125        elif icod == -1:
126            resList.append((None,iref))
127        elif icod == -2:
128            break
129    return sInt,resList
130
131def ComputeFobsSqCW(refl,iref):
132    yp = np.zeros(len(x)) # not masked
133    sInt = 0
134    refl8im = 0
135    Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
136    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
137    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
138    iFin2 = iFin
139    if not iBeg+iFin:       #peak below low limit - skip peak
140        return 0
141    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
142        return -1
143    elif not iBeg-iFin:     #peak above high limit - done
144        return -2
145    elif iBeg < iFin:
146        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(
147            refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin])
148        sInt = refl[11+im]*refl[9+im]
149        if Ka2:
150            pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
151            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
152            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
153            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
154            if iFin2 > iBeg2: 
155                yp[iBeg2:iFin2] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(
156                    pos2,refl[6+im],refl[7+im],shl,x[iBeg2:iFin2])
157                sInt *= 1.+kRatio
158    refl8im = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
159    return refl8im,sInt
160
161def ComputeFobsSqTOF(refl,iref):
162    yp = np.zeros(len(x)) # not masked
163    refl8im = 0
164    Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
165    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
166    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
167    if not iBeg+iFin:       #peak below low limit - skip peak
168        return 0
169    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
170        return -1
171    elif not iBeg-iFin:     #peak above high limit - done
172        return -2
173    if iBeg < iFin:
174        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(
175            refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])
176    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
177    return refl8im,refl[11+im]*refl[9+im]
178
179def ComputeFobsSqPink(refl,iref):
180    yp = np.zeros(len(x)) # not masked
181    refl8im = 0
182    Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.)
183    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
184    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
185    if not iBeg+iFin:       #peak below low limit - skip peak
186        return 0
187    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
188        return -1
189    elif not iBeg-iFin:     #peak above high limit - done
190        return -2
191    if iBeg < iFin:
192        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(
193            refl[5+im],refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.,x[iBeg:iFin])
194    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
195    return refl8im,refl[11+im]*refl[9+im]
196
197################################################################################
198# Powder Profile computation
199################################################################################       
200def InitPwdrProfGlobals(im1,shl1,x1):
201    '''Initialize for the computation of Fobs Squared for powder histograms.
202    Puts lots of junk into the global namespace in this module.
203    '''
204    global im,shl,x
205    im = im1
206    shl = shl1
207    x = ma.getdata(x1)
208    global cw
209    cw = np.diff(x)
210    cw = np.append(cw,cw[-1])
211    # create local copies of ycalc array
212    global yc
213    yc = np.zeros_like(x)
214
215
216def ComputePwdrProfCW(profList):
217    'Compute the peaks profile for a set of CW peaks and add into the yc array'
218    for pos,refl,iBeg,iFin,kRatio in profList:
219        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(
220            pos,refl[6+im],refl[7+im],shl,x[iBeg:iFin])
221    return yc
222
223def ComputePwdrProfTOF(profList):
224    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
225    for pos,refl,iBeg,iFin in profList:
226        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(
227            pos,refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])/cw[iBeg:iFin]
228    return yc
229   
230def ComputePwdrProfPink(profList):
231    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
232    for pos,refl,iBeg,iFin in profList:
233        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(
234            pos,refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.,x[iBeg:iFin])/cw[iBeg:iFin]
235    return yc
Note: See TracBrowser for help on using the repository browser.