source: trunk/GSASIImpsubs.py @ 4521

Last change on this file since 4521 was 4521, checked in by vondreele, 15 months 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.