Ignore:
Timestamp:
Jul 25, 2017 6:21:06 PM (4 years ago)
Author:
toby
Message:

Finally: multiprocessing for refinement

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branch/2frame/GSASIImpsubs.py

    r2941 r2944  
    2828import numpy as np
    2929import numpy.ma as ma
     30import numpy.linalg as nl
    3031import GSASIIpath
    3132GSASIIpath.SetVersionNumber("$Revision: 2895 $")
     
    3334import GSASIIstrMath as G2stMth
    3435
     36sind = lambda x: np.sin(x*np.pi/180.)
    3537cosd = lambda x: np.cos(x*np.pi/180.)
    3638tand = lambda x: np.tan(x*np.pi/180.)
     39#asind = lambda x: 180.*np.arcsin(x)/np.pi
    3740#acosd = lambda x: 180.*np.arccos(x)/np.pi
    3841#atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
     
    4043ncores = None
    4144
    42 def InitMP():
     45def InitMP(allowMP=True):
    4346    '''Called in routines to initialize use of Multiprocessing
    4447    '''
     
    4649    if ncores is not None: return
    4750    useMP = False
     51    if not allowMP:
     52        print('Multiprocessing disabled')
     53        ncores = 0
     54        return
    4855    ncores = GSASIIpath.GetConfigValue('Multiprocessing_cores',-1)
    4956    if ncores < 0: ncores = mp.cpu_count()
     
    5966def InitDerivGlobals(im1,calcControls1,SGMT1,hfx1,phfx1,pfx1,G1,GB1,g1,SGData1,
    6067                     parmDict1,wave1,shl1,x1,cw1,Ka21,A1,varylist1,dependentVars1,
    61                      dFdvDict1,lamRatio1,kRatio1,doPawley1):
     68                     dFdvDict1,lamRatio1,kRatio1,doPawley1,pawleyLookup1):
    6269    '''Initialize for the computation of derivatives. Puts lots of junk into the global
    6370    namespace in this module, including the arrays for derivatives (when needed.)
    6471    '''
    6572    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
     73    global varylist,dependentVars,dFdvDict,lamRatio,kRatio,doPawley,pawleyLookup
    6774    im = im1
    6875    calcControls = calcControls1
     
    8895    kRatio = kRatio1
    8996    doPawley = doPawley1
     97    pawleyLookup = pawleyLookup1
    9098    # determine the parameters that will have derivatives computed only at end
    9199    global nonatomvarylist
     
    164172    dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = G2stMth.GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    165173    pos = refl[5+im]
     174    calcKa2 = False
    166175    if 'C' in calcControls[hfx+'histType']:
    167176        tanth = tand(pos/2.0)
    168177        costh = cosd(pos/2.0)
    169         calcKa2 = False
    170178        if Ka2:
    171179            pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
     
    222230            names.update({hfx+'Absorption':[dFdAb,'int'],})
    223231    else:   #'T'OF
    224         dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
     232        dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = G2stMth.GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
    225233        names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
    226234            hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
     
    345353           #if calcKa2:
    346354           #    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr2
     355
     356################################################################################
     357# Fobs Squared computation
     358################################################################################       
     359#x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2
     360def InitFobsSqGlobals(x1,ratio1,shl1,xB1,xF1,im1,lamRatio1,kRatio1,xMask1,Ka21):
     361    '''Initialize for the computation of Fobs Squared for powder histograms.
     362    Puts lots of junk into the global namespace in this module.
     363    '''
     364    global x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2
     365    x = ma.getdata(x1)
     366    ratio = ratio1
     367    shl = shl1
     368    xB = xB1
     369    xF = xF1
     370    im = im1
     371    lamRatio = lamRatio1
     372    kRatio = kRatio1
     373    xMask = xMask1
     374    Ka2 = Ka21
     375
     376def ComputeFobsSqCWbatch(profList):
     377    sInt = 0
     378    resList = []
     379    for refl,iref in profList:
     380        icod = ComputeFobsSqCW(refl,iref)
     381        if type(icod) is tuple:
     382            resList.append((icod[0],iref))
     383            sInt += icod[1]
     384        elif icod == -1:
     385            res.append((None,iref))
     386        elif icod == -2:
     387            break
     388    return sInt,resList
     389
     390def ComputeFobsSqTOFbatch(profList):
     391    sInt = 0
     392    resList = []
     393    for refl,iref in profList:
     394        icod = ComputeFobsSqTOF(refl,iref)
     395        if type(icod) is tuple:
     396            resList.append((icod[0],iref))
     397            sInt += icod[1]
     398        elif icod == -1:
     399            res.append((None,iref))
     400        elif icod == -2:
     401            break
     402    return sInt,resList
     403       
     404def ComputeFobsSqCW(refl,iref):
     405    yp = np.zeros(len(x)) # not masked
     406    sInt = 0
     407    refl8im = 0
     408    Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
     409    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
     410    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
     411    iFin2 = iFin
     412    if not iBeg+iFin:       #peak below low limit - skip peak
     413        return 0
     414    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
     415        return -1
     416    elif not iBeg-iFin:     #peak above high limit - done
     417        return -2
     418    elif iBeg < iFin:
     419        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(
     420            refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin])
     421        sInt = refl[11+im]*refl[9+im]
     422        if Ka2:
     423            pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     424            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
     425            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
     426            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
     427            if iFin2 > iBeg2:
     428                yp[iBeg2:iFin2] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(
     429                    pos2,refl[6+im],refl[7+im],shl,x[iBeg2:iFin2])
     430                sInt *= 1.+kRatio
     431    refl8im = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
     432    return refl8im,sInt
     433
     434def ComputeFobsSqTOF(refl,iref):
     435    yp = np.zeros(len(x)) # not masked
     436    refl8im = 0
     437    Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
     438    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
     439    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
     440    if not iBeg+iFin:       #peak below low limit - skip peak
     441        return 0
     442    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
     443        return -1
     444    elif not iBeg-iFin:     #peak above high limit - done
     445        return -2
     446    if iBeg < iFin:
     447        yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(
     448            refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])
     449    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
     450    return refl8im,refl[11+im]*refl[9+im]
     451################################################################################
     452# Powder Profile computation
     453################################################################################       
     454def InitPwdrProfGlobals(im1,shl1,x1):
     455    '''Initialize for the computation of Fobs Squared for powder histograms.
     456    Puts lots of junk into the global namespace in this module.
     457    '''
     458    global im,shl,x
     459    im = im1
     460    shl = shl1
     461    x = ma.getdata(x1)
     462    global cw
     463    cw = np.diff(x)
     464    cw = np.append(cw,cw[-1])
     465    # create local copies of ycalc array
     466    if useMP:
     467        global yc
     468        yc = np.zeros_like(x1)
     469
     470
     471def ComputePwdrProfCW(profList):
     472    'Compute the peaks profile for a set of CW peaks and add into the yc array'
     473    for pos,refl,iBeg,iFin,kRatio in profList:
     474        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(
     475            pos,refl[6+im],refl[7+im],shl,x[iBeg:iFin])
     476    return yc
     477
     478def ComputePwdrProfTOF(profList):
     479    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
     480    for pos,refl,iBeg,iFin in profList:
     481        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(
     482            pos,refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])/cw[iBeg:iFin]
     483    return yc
     484   
Note: See TracChangeset for help on using the changeset viewer.