source: trunk/GSASIImpsubs.py @ 3866

Last change on this file since 3866 was 3866, checked in by vondreele, 3 years ago

fix problem with multicore range error - need integer divide in ResetMP

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 6.8 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: 2019-04-01 19:13:10 +0000 (Mon, 01 Apr 2019) $
22# $Author: vondreele $
23# $Revision: 3866 $
24# $URL: trunk/GSASIImpsubs.py $
25# $Id: GSASIImpsubs.py 3866 2019-04-01 19:13:10Z 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: 3866 $")
33import GSASIIpwd as G2pwd
34
35sind = lambda x: np.sin(x*np.pi/180.)
36cosd = lambda x: np.cos(x*np.pi/180.)
37tand = lambda x: np.tan(x*np.pi/180.)
38#asind = lambda x: 180.*np.arcsin(x)/np.pi
39#acosd = lambda x: 180.*np.arccos(x)/np.pi
40#atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
41   
42ncores = None
43
44def ResetMP():
45    '''Call after changing Config var 'Multiprocessing_cores' to force a resetting
46    of the useMP from the parameter.
47    '''
48    global ncores
49    ncores = None
50   
51def InitMP(allowMP=True):
52    '''Called to initialize use of Multiprocessing
53    '''
54    global useMP,ncores
55    if ncores is not None: return useMP,ncores
56    useMP = False
57    if not allowMP:
58        print('Multiprocessing disabled')
59        ncores = 0
60        return useMP,ncores
61    ncores = GSASIIpath.GetConfigValue('Multiprocessing_cores',0)
62    if ncores < 0: ncores = mp.cpu_count()//2
63    if ncores > 1:
64        useMP = True
65    if useMP:
66        print('Multiprocessing with {} cores enabled'.format(ncores))
67    return useMP,ncores
68
69################################################################################
70# Fobs Squared computation
71################################################################################       
72def InitFobsSqGlobals(x1,ratio1,shl1,xB1,xF1,im1,lamRatio1,kRatio1,xMask1,Ka21):
73    '''Initialize for the computation of Fobs Squared for powder histograms.
74    Puts lots of junk into the global namespace in this module.
75    '''
76    global x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2
77    x = ma.getdata(x1)
78    ratio = ratio1
79    shl = shl1
80    xB = xB1
81    xF = xF1
82    im = im1
83    lamRatio = lamRatio1
84    kRatio = kRatio1
85    xMask = xMask1
86    Ka2 = Ka21
87
88def ComputeFobsSqCWbatch(profList):
89    sInt = 0
90    resList = []
91    for refl,iref in profList:
92        icod = ComputeFobsSqCW(refl,iref)
93        if type(icod) is tuple:
94            resList.append((icod[0],iref))
95            sInt += icod[1]
96        elif icod == -1:
97            resList.append((None,iref))
98        elif icod == -2:
99            break
100    return sInt,resList
101
102def ComputeFobsSqTOFbatch(profList):
103    sInt = 0
104    resList = []
105    for refl,iref in profList:
106        icod = ComputeFobsSqTOF(refl,iref)
107        if type(icod) is tuple:
108            resList.append((icod[0],iref))
109            sInt += icod[1]
110        elif icod == -1:
111            resList.append((None,iref))
112        elif icod == -2:
113            break
114    return sInt,resList
115       
116def ComputeFobsSqCW(refl,iref):
117    yp = np.zeros(len(x)) # not masked
118    sInt = 0
119    refl8im = 0
120    Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
121    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
122    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
123    iFin2 = iFin
124    if not iBeg+iFin:       #peak below low limit - skip peak
125        return 0
126    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
127        return -1
128    elif not iBeg-iFin:     #peak above high limit - done
129        return -2
130    elif iBeg < iFin:
131        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(
132            refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin])
133        sInt = refl[11+im]*refl[9+im]
134        if Ka2:
135            pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
136            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
137            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
138            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
139            if iFin2 > iBeg2: 
140                yp[iBeg2:iFin2] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(
141                    pos2,refl[6+im],refl[7+im],shl,x[iBeg2:iFin2])
142                sInt *= 1.+kRatio
143    refl8im = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
144    return refl8im,sInt
145
146def ComputeFobsSqTOF(refl,iref):
147    yp = np.zeros(len(x)) # not masked
148    refl8im = 0
149    Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
150    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
151    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
152    if not iBeg+iFin:       #peak below low limit - skip peak
153        return 0
154    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
155        return -1
156    elif not iBeg-iFin:     #peak above high limit - done
157        return -2
158    if iBeg < iFin:
159        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(
160            refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])
161    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
162    return refl8im,refl[11+im]*refl[9+im]
163################################################################################
164# Powder Profile computation
165################################################################################       
166def InitPwdrProfGlobals(im1,shl1,x1):
167    '''Initialize for the computation of Fobs Squared for powder histograms.
168    Puts lots of junk into the global namespace in this module.
169    '''
170    global im,shl,x
171    im = im1
172    shl = shl1
173    x = ma.getdata(x1)
174    global cw
175    cw = np.diff(x)
176    cw = np.append(cw,cw[-1])
177    # create local copies of ycalc array
178    global yc
179    yc = np.zeros_like(x)
180
181
182def ComputePwdrProfCW(profList):
183    'Compute the peaks profile for a set of CW peaks and add into the yc array'
184    for pos,refl,iBeg,iFin,kRatio in profList:
185        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(
186            pos,refl[6+im],refl[7+im],shl,x[iBeg:iFin])
187    return yc
188
189def ComputePwdrProfTOF(profList):
190    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
191    for pos,refl,iBeg,iFin in profList:
192        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(
193            pos,refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])/cw[iBeg:iFin]
194    return yc
195   
Note: See TracBrowser for help on using the repository browser.