Changeset 2948


Ignore:
Timestamp:
Jul 29, 2017 8:06:56 PM (6 years ago)
Author:
toby
Message:

return to version [2935] of GSASIIstrMath.py w/o MP

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branch/2frame/GSASIIstrMath.py

    r2947 r2948  
    1313import time
    1414import copy
    15 import multiprocessing as mp
    1615import numpy as np
    1716import numpy.ma as ma
    … …  
    2726import GSASIImath as G2mth
    2827import GSASIIobj as G2obj
    29 import GSASIImpsubs as G2mp
    3028
    3129sind = lambda x: np.sin(x*np.pi/180.)
    … …  
    723721#reflection processing begins here - big arrays!
    724722    iBeg = 0
    725 #    time0 = time.time()
     723    time0 = time.time()
    726724    while iBeg < nRef:
    727725        iFin = min(iBeg+blkSize,nRef)
    … …  
    862860    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    863861        Flack = 1.-2.*parmDict[phfx+'Flack']
    864 #    time0 = time.time()
     862    time0 = time.time()
    865863#reflection processing begins here - big arrays!
    866864    iBeg = 0
    … …  
    29262924def GetFobsSq(Histograms,Phases,parmDict,calcControls):
    29272925    'Compute the observed structure factors for Powder histograms'
    2928     #if GSASIIpath.GetConfigValue('debug'):
    2929     #    starttime = time.time(); print 'start GetFobsSq'
    2930     G2mp.InitMP()
     2926    #starttime = time.time(); print 'start GetFobsSq'
    29312927    histoList = Histograms.keys()
    29322928    histoList.sort()
    2933     Ka2 = shl = lamRatio = kRatio = None
    29342929    for histogram in histoList:
    29352930        if 'PWDR' in histogram[:4]:
    … …  
    29532948            ymb = np.where(ymb,ymb,1.0)
    29542949            ycmb = np.array(yc-yb)
    2955             ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
     2950            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
    29562951            refLists = Histogram['Reflection Lists']
    29572952            for phase in refLists:
    … …  
    29652960                phfx = '%d:%d:'%(pId,hId)
    29662961                refDict = refLists[phase]
    2967                 sumInt = 0.0
    2968                 nExcl = 0
    2969                 useMP = G2mp.useMP and len(refDict['RefList']) > 100
    2970                 if useMP: # yes, create a set of initialized Python processes
    2971                     MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
    2972                                     [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
    2973                     profArgs = [[] for i in range(G2mp.ncores)]
    2974                 else:
    2975                     G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
    2976                 if 'C' in calcControls[hfx+'histType']:
    2977                     # are we multiprocessing?
    2978                     for iref,refl in enumerate(refDict['RefList']):
    2979                         if useMP:
    2980                             profArgs[iref%G2mp.ncores].append((refl,iref))
    2981                         else:
    2982                             icod= G2mp.ComputeFobsSqCW(refl,iref)
    2983                             if type(icod) is tuple:
    2984                                 refl[8+im] = icod[0]
    2985                                 sumInt += icod[1]
    2986                                 if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
    2987                             elif icod == -1:
    2988                                 refl[3+im] *= -1
    2989                                 nExcl += 1
    2990                             elif icod == -2:
    2991                                 break
    2992                     if useMP:
    2993                         for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
    2994                             sumInt += sInt
    2995                             for refl8im,irefl in resList:
    2996                                 if refl8im is None:
    2997                                     refDict['RefList'][irefl][3+im] *= -1
    2998                                     nExcl += 1
    2999                                 else:
    3000                                     refDict['RefList'][irefl][8+im] = refl8im
    3001                                     if parmDict[phfx+'LeBail']:
    3002                                         refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
    3003                         MPpool.terminate()
    3004                 elif 'T' in calcControls[hfx+'histType']:
    3005                     for iref,refl in enumerate(refDict['RefList']):
    3006                         if useMP:
    3007                             profArgs[iref%G2mp.ncores].append((refl,iref))
    3008                         else:
    3009                             icod= G2mp.ComputeFobsSqTOF(refl,iref)
    3010                             if type(icod) is tuple:
    3011                                 refl[8+im] = icod[0]
    3012                                 sumInt += icod[1]
    3013                                 if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
    3014                             elif icod == -1:
    3015                                 refl[3+im] *= -1
    3016                                 nExcl += 1
    3017                             elif icod == -2:
    3018                                 break
    3019                     if useMP:
    3020                         for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
    3021                             sumInt += sInt
    3022                             for refl8im,irefl in resList:
    3023                                 if refl8im is None:
    3024                                     refDict['RefList'][irefl][3+im] *= -1
    3025                                     nExcl += 1
    3026                                 else:
    3027                                     refDict['RefList'][irefl][8+im] = refl8im
    3028                                     if parmDict[phfx+'LeBail']:
    3029                                         refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
    30302962                sumFo = 0.0
    30312963                sumdF = 0.0
    30322964                sumFosq = 0.0
    30332965                sumdFsq = 0.0
    3034                 for iref,refl in enumerate(refDict['RefList']):
     2966                sumInt = 0.0
     2967                nExcl = 0
     2968                for refl in refDict['RefList']:
     2969                    if 'C' in calcControls[hfx+'histType']:
     2970                        yp = np.zeros_like(yb)
     2971                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
     2972                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
     2973                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
     2974                        iFin2 = iFin
     2975                        if not iBeg+iFin:       #peak below low limit - skip peak
     2976                            continue
     2977                        if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
     2978                            refl[3+im] *= -1
     2979                            nExcl += 1
     2980                            continue
     2981                        elif not iBeg-iFin:     #peak above high limit - done
     2982                            break
     2983                        elif iBeg < iFin:
     2984                            yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
     2985                            sumInt += refl[11+im]*refl[9+im]
     2986                            if Ka2:
     2987                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     2988                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
     2989                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
     2990                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
     2991                                if iFin2 > iBeg2:
     2992                                    yp[iBeg2:iFin2] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))        #and here
     2993                                    sumInt += refl[11+im]*refl[9+im]*kRatio
     2994                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
     2995                            if parmDict[phfx+'LeBail']:
     2996                                refl[9+im] = refl[8+im]
     2997                               
     2998                    elif 'T' in calcControls[hfx+'histType']:
     2999                        yp = np.zeros_like(yb)
     3000                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
     3001                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
     3002                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
     3003                        if not iBeg+iFin:       #peak below low limit - skip peak
     3004                            continue
     3005                        if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
     3006                            refl[3+im] *= -1
     3007                            nExcl += 1
     3008                            continue
     3009                        elif not iBeg-iFin:     #peak above high limit - done
     3010                            break
     3011                        if iBeg < iFin:
     3012                            yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))  #>90% of time spent here
     3013                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
     3014                            if parmDict[phfx+'LeBail']:
     3015                                refl[9+im] = refl[8+im]
     3016                            sumInt += refl[11+im]*refl[9+im]
    30353017                    Fo = np.sqrt(np.abs(refl[8+im]))
    30363018                    Fc = np.sqrt(np.abs(refl[9]+im))
    … …  
    30513033            Histogram = Histograms[histogram]
    30523034            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
    3053     #if GSASIIpath.GetConfigValue('debug'):
    3054     #    print 'end GetFobsSq t=',time.time()-starttime
    3055     #    print "Rf^2",[(i,Histogram['Residuals'][i]) for i in Histogram['Residuals'] if 'Rf^2' in i]
    3056 
     3035    #print 'end GetFobsSq t=',time.time()-starttime
     3036               
    30573037def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
    30583038    'Computes the powder pattern for a histogram based on contributions from all used phases'
    3059     #if GSASIIpath.GetConfigValue('debug'):
    3060     #    starttime = time.time(); print 'start getPowderProfile'
    3061     G2mp.InitMP()   
    3062            
     3039   
     3040    #starttime = time.time(); print 'start getPowderProfile'
     3041   
    30633042    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
    30643043        U = parmDict[hfx+'U']
    … …  
    30883067        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
    30893068        return alp,bet
    3090    
     3069       
    30913070    hId = Histogram['hId']
    30923071    hfx = ':%d:'%(hId)
    … …  
    30963075    cw = np.diff(x)
    30973076    cw = np.append(cw,cw[-1])
    3098    
    3099     shl = wave = kRatio = lamRatio = None # variables that are passed and must be initialized even when not needed
     3077       
    31003078    if 'C' in calcControls[hfx+'histType']:   
    31013079        shl = max(parmDict[hfx+'SH/L'],0.002)
    … …  
    31373115                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    31383116        badPeak = False
    3139 
    3140         # are we multiprocessing?
    3141         useMP = G2mp.useMP and len(refDict['RefList']) > 100
    3142         if useMP: # yes, create a set of initialized Python processes
    3143             MPpool = mp.Pool(G2mp.ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
    3144             profArgs = [[] for i in range(G2mp.ncores)]
    3145         if 'C' in calcControls[hfx+'histType']:
    3146             for iref,refl in enumerate(refDict['RefList']):
     3117        for iref,refl in enumerate(refDict['RefList']):
     3118            if 'C' in calcControls[hfx+'histType']:
    31473119                if im:
    31483120                    h,k,l,m = refl[:4]
    … …  
    31763148                    badPeak = True
    31773149                    continue
    3178                 if useMP:
    3179                     profArgs[iref%G2mp.ncores].append((refl[5+im],refl,iBeg,iFin,1.))
    3180                 else:
    3181                     yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(
    3182                         refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
     3150                yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
    31833151                if Ka2:
    31843152                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    … …  
    31923160                    elif iBeg > iFin:   #bad peak coeff - skip
    31933161                        continue
    3194                     if useMP:
    3195                         profArgs[iref%G2mp.ncores].append((pos2,refl,iBeg,iFin,kRatio))
    3196                     else:
    3197                         yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(
    3198                             pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
    3199         elif 'T' in calcControls[hfx+'histType']:
    3200             for iref,refl in enumerate(refDict['RefList']):
     3162                    yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))        #and here
     3163            elif 'T' in calcControls[hfx+'histType']:
    32013164                h,k,l = refl[:3]
    32023165                Uniq = np.inner(refl[:3],SGMT)
    3203                 refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
     3166                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position - #TODO - what about tabluated offset?
    32043167                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
    32053168#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
    32063169                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
    3207                 refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
     3170                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)             #TODO - skip if alp, bet tabulated?
    32083171                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    32093172                refl[11+im] *= Vst*Lorenz
    … …  
    32283191                    badPeak = True
    32293192                    continue
    3230                 if useMP:
    3231                     profArgs[iref%G2mp.ncores].append((refl[5+im],refl,iBeg,iFin))
    3232                 else:
    3233                     yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
    3234         if useMP and 'C' in calcControls[hfx+'histType']:
    3235             for y in MPpool.imap_unordered(G2mp.ComputePwdrProfCW,profArgs):
    3236                 yc += y
    3237             MPpool.terminate()
    3238         elif useMP:
    3239             for y in MPpool.imap_unordered(G2mp.ComputePwdrProfTOF,profArgs):
    3240                 yc += y
    3241             MPpool.terminate()
     3193                yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
     3194#        print 'profile calc time: %.3fs'%(time.time()-time0)
    32423195    if badPeak:
    32433196        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped'
    3244     #if GSASIIpath.GetConfigValue('debug'):
    3245     #    print 'end getPowderProfile t=',time.time()-starttime
     3197    #print 'end getPowderProfile t=',time.time()-starttime
    32463198    return yc,yb
    32473199   
    … …  
    32533205    #    starttime = time.time()
    32543206    #    print 'starting getPowderProfileDerv'
    3255     G2mp.InitMP()   
     3207   
     3208    def cellVaryDerv(pfx,SGData,dpdA):
     3209        if SGData['SGLaue'] in ['-1',]:
     3210            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
     3211                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
     3212        elif SGData['SGLaue'] in ['2/m',]:
     3213            if SGData['SGUniq'] == 'a':
     3214                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
     3215            elif SGData['SGUniq'] == 'b':
     3216                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
     3217            else:
     3218                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
     3219        elif SGData['SGLaue'] in ['mmm',]:
     3220            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
     3221        elif SGData['SGLaue'] in ['4/m','4/mmm']:
     3222            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
     3223        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
     3224            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
     3225        elif SGData['SGLaue'] in ['3R', '3mR']:
     3226            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
     3227        elif SGData['SGLaue'] in ['m3m','m3']:
     3228            return [[pfx+'A0',dpdA[0]]]
     3229           
    32563230    # create a list of dependent variables and set up a dictionary to hold their derivatives
    32573231    dependentVars = G2mv.GetDependentVars()
    … …  
    32883262    cw = np.append(cw,cw[-1])
    32893263    Ka2 = False #also for TOF!
    3290     dFdvDict = shl = wave = kRatio = lamRatio = None # variables that are passed and must be initialized even when not needed
    32913264    if 'C' in calcControls[hfx+'histType']:   
    32923265        shl = max(parmDict[hfx+'SH/L'],0.002)
    … …  
    32983271        else:
    32993272            wave = parmDict[hfx+'Lam']
     3273    #print '#1 getPowderProfileDerv t=',time.time()-starttime
    33003274    for phase in Histogram['Reflection Lists']:
    33013275        refDict = Histogram['Reflection Lists'][phase]
    … …  
    33273301#            print 'sf-derv time %.3fs'%(time.time()-time0)
    33283302            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
    3329 
    3330         # are we multiprocessing?
    3331         useMP = G2mp.useMP and len(refDict['RefList']) > 100
    3332         if useMP: # yes, create a set of initialized Python processes
    3333             MPpool = mp.Pool(G2mp.ncores,G2mp.InitDerivGlobals,
    3334                             [im,calcControls,SGMT,hfx,phfx,pfx,G,GB,g,SGData,parmDict,
    3335                                wave,shl,x,cw,Ka2,A,varylist,dependentVars,dFdvDict,
    3336                                lamRatio,kRatio,Phase['General'].get('doPawley'),
    3337                                pawleyLookup
    3338                             ])
    3339             derivArgs = [[] for i in range(G2mp.ncores)]
    3340         else: # no, initialize the G2mp module
    3341             G2mp.InitDerivGlobals(im,calcControls,SGMT,hfx,phfx,pfx,G,GB,g,SGData,parmDict,
    3342                               wave,shl,x,cw,Ka2,A,varylist,dependentVars,dFdvDict,
    3343                               lamRatio,kRatio,Phase['General'].get('doPawley'),
    3344                               pawleyLookup)
    3345         # loop over all reflections
     3303        #print '#2 getPowderProfileDerv t=',time.time()-starttime
     3304        # determine the parameters that will have derivatives computed only at end
     3305        nonatomvarylist = []
     3306        for name in varylist:
     3307            if '::RBV;' not in name:
     3308                try:
     3309                    aname = name.split(pfx)[1][:2]
     3310                    if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
     3311                        'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
     3312                except IndexError:
     3313                    continue
     3314            nonatomvarylist.append(name)
     3315        nonatomdependentVars = []
     3316        for name in dependentVars:
     3317            if '::RBV;' not in name:
     3318                try:
     3319                    aname = name.split(pfx)[1][:2]
     3320                    if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
     3321                        'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
     3322                except IndexError:
     3323                    continue
     3324            nonatomdependentVars.append(name)
     3325        #timelist =  10*[0.0]
     3326        #timestart = 10*[0.0]
     3327        #==========================================================================================
     3328        #==========================================================================================
    33463329        for iref,refl in enumerate(refDict['RefList']):
     3330            #timestart[0] = time.time()
     3331            if im:
     3332                h,k,l,m = refl[:4]
     3333            else:
     3334                h,k,l = refl[:3]
     3335            Uniq = np.inner(refl[:3],SGMT)
     3336            if 'T' in calcControls[hfx+'histType']:
     3337                wave = refl[14+im]
     3338            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    33473339            if 'C' in calcControls[hfx+'histType']:        #CW powder
    33483340                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
    … …  
    33553347            elif not iBeg-iFin:     #peak above high limit - done
    33563348                break
    3357             if useMP:
    3358                 derivArgs[iref%G2mp.ncores].append((refl,iref,fmin,fmax,iBeg,iFin))
    3359             else:
    3360                 result = G2mp.ComputeDeriv(refl,iref,fmin,fmax,iBeg,iFin,dMdv,depDerivDict)
    3361                 if result: break
    3362         if useMP:
    3363             for dm,depd in MPpool.imap_unordered(G2mp.ComputeDerivMPbatch,derivArgs):
    3364                 dMdv += dm
    3365                 for i in depDerivDict: depDerivDict[i] += depd[i]
    3366             MPpool.terminate()
     3349            pos = refl[5+im]
     3350            #itim=0;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3351            if 'C' in calcControls[hfx+'histType']:
     3352                tanth = tand(pos/2.0)
     3353                costh = cosd(pos/2.0)
     3354                lenBF = iFin-iBeg
     3355                dMdpk = np.zeros(shape=(6,lenBF))
     3356                dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
     3357                for i in range(5):
     3358                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
     3359                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
     3360                if Ka2:
     3361                    pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
     3362                    iBeg2 = np.searchsorted(x,pos2-fmin)
     3363                    iFin2 = np.searchsorted(x,pos2+fmax)
     3364                    if iBeg2-iFin2:
     3365                        lenBF2 = iFin2-iBeg2
     3366                        dMdpk2 = np.zeros(shape=(6,lenBF2))
     3367                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
     3368                        for i in range(5):
     3369                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
     3370                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
     3371                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
     3372            else:   #'T'OF
     3373                lenBF = iFin-iBeg
     3374                if lenBF < 0:   #bad peak coeff
     3375                    break
     3376                dMdpk = np.zeros(shape=(6,lenBF))
     3377                dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
     3378                for i in range(6):
     3379                    dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
     3380                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
     3381            #itim=1;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3382            if Phase['General'].get('doPawley'):
     3383                dMdpw = np.zeros(len(x))
     3384                try:
     3385                    if im:
     3386                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
     3387                    else:
     3388                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
     3389                    idx = varylist.index(pIdx)
     3390                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
     3391                    if Ka2: #not for TOF either
     3392                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
     3393                    dMdv[idx] = dMdpw
     3394                except: # ValueError:
     3395                    pass
     3396            if 'C' in calcControls[hfx+'histType']:
     3397                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
     3398                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
     3399                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
     3400                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
     3401                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
     3402                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
     3403                    hfx+'DisplaceY':[dpdY,'pos'],}
     3404                if 'Bragg' in calcControls[hfx+'instType']:
     3405                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
     3406                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
     3407                else:
     3408                    names.update({hfx+'Absorption':[dFdAb,'int'],})
     3409            else:   #'T'OF
     3410                dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
     3411                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
     3412                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
     3413                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
     3414                    hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
     3415                    hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
     3416                    hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
     3417                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
     3418            #itim=2;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3419            for name in names:
     3420                item = names[name]
     3421                if name in varylist:
     3422                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
     3423                    if Ka2 and iFin2-iBeg2:
     3424                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
     3425                elif name in dependentVars:
     3426                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
     3427                    if Ka2 and iFin2-iBeg2:
     3428                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
     3429            for iPO in dIdPO:
     3430                if iPO in varylist:
     3431                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
     3432                    if Ka2 and iFin2-iBeg2:
     3433                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
     3434                elif iPO in dependentVars:
     3435                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
     3436                    if Ka2 and iFin2-iBeg2:
     3437                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
     3438            #itim=3;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3439            for i,name in enumerate(['omega','chi','phi']):
     3440                aname = pfx+'SH '+name
     3441                if aname in varylist:
     3442                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
     3443                    if Ka2 and iFin2-iBeg2:
     3444                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
     3445                elif aname in dependentVars:
     3446                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
     3447                    if Ka2 and iFin2-iBeg2:
     3448                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
     3449            for iSH in dFdODF:
     3450                if iSH in varylist:
     3451                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
     3452                    if Ka2 and iFin2-iBeg2:
     3453                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
     3454                elif iSH in dependentVars:
     3455                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
     3456                    if Ka2 and iFin2-iBeg2:
     3457                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
     3458            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
     3459            for name,dpdA in cellDervNames:
     3460                if name in varylist:
     3461                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
     3462                    if Ka2 and iFin2-iBeg2:
     3463                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
     3464                elif name in dependentVars: #need to scale for mixed phase constraints?
     3465                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
     3466                    if Ka2 and iFin2-iBeg2:
     3467                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
     3468            dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
     3469            for name in dDijDict:
     3470                if name in varylist:
     3471                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
     3472                    if Ka2 and iFin2-iBeg2:
     3473                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
     3474                elif name in dependentVars:
     3475                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
     3476                    if Ka2 and iFin2-iBeg2:
     3477                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
     3478            #itim=4;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3479            for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
     3480                if name in varylist:
     3481                    dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
     3482                    if Ka2 and iFin2-iBeg2:
     3483                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
     3484                elif name in dependentVars:
     3485                    depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
     3486                    if Ka2 and iFin2-iBeg2:
     3487                        depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
     3488            if 'C' in calcControls[hfx+'histType']:
     3489                sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
     3490            else:   #'T'OF
     3491                sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
     3492            for name in gamDict:
     3493                if name in varylist:
     3494                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
     3495                    if Ka2 and iFin2-iBeg2:
     3496                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
     3497                elif name in dependentVars:
     3498                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
     3499                    if Ka2 and iFin2-iBeg2:
     3500                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
     3501            for name in sigDict:
     3502                if name in varylist:
     3503                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
     3504                    if Ka2 and iFin2-iBeg2:
     3505                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
     3506                elif name in dependentVars:
     3507                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
     3508                    if Ka2 and iFin2-iBeg2:
     3509                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
     3510            for name in ['BabA','BabU']:
     3511                if refl[9+im]:
     3512                    if phfx+name in varylist:
     3513                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
     3514                        if Ka2 and iFin2-iBeg2:
     3515                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]
     3516                    elif phfx+name in dependentVars:                   
     3517                        depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
     3518                        if Ka2 and iFin2-iBeg2:
     3519                            depDerivDict[phfx+name][iBeg2:iFin2] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]                 
     3520            #itim=5;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3521            if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
     3522                #do atom derivatives -  for RB,F,X & U so far - how do I scale mixed phase constraints?
     3523                corr = 0.
     3524                corr2 = 0.
     3525                if refl[9+im]:             
     3526                    corr = dervDict['int']/refl[9+im]
     3527                    #if Ka2 and iFin2-iBeg2:
     3528                    #    corr2 = dervDict2['int']/refl[9+im]
     3529                #itim=6;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3530                for name in nonatomvarylist:
     3531                    dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
     3532                    if Ka2 and iFin2-iBeg2:
     3533                       dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     3534                #itim=7;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3535                for name in nonatomdependentVars:
     3536                   depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
     3537                   if Ka2 and iFin2-iBeg2:
     3538                       depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     3539                #itim=8;timelist[itim] += time.time()-timestart[itim]   
     3540    #        print 'profile derv time: %.3fs'%(time.time()-time0)
    33673541    # now process derivatives in constraints
     3542    #print '#3 getPowderProfileDerv t=',time.time()-starttime
     3543    #print timelist,sum(timelist)
    33683544    dMdv = ma.array(dMdv,mask=np.outer(np.ones(len(varylist)),ma.getmaskarray(x)))      #x is a MaskedArray!
    33693545    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
    … …  
    38454021        M = np.concatenate((M,np.sqrt(pWt)*pVals))
    38464022    return M
    3847 
     4023                       
Note: See TracChangeset for help on using the changeset viewer.