Ignore:
Timestamp:
Jul 29, 2017 8:06:56 PM (4 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.