 Timestamp:
 Jul 29, 2017 8:06:56 PM (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branch/2frame/GSASIIstrMath.py
r2947 r2948 13 13 import time 14 14 import copy 15 import multiprocessing as mp16 15 import numpy as np 17 16 import numpy.ma as ma … … 27 26 import GSASIImath as G2mth 28 27 import GSASIIobj as G2obj 29 import GSASIImpsubs as G2mp30 28 31 29 sind = lambda x: np.sin(x*np.pi/180.) … … 723 721 #reflection processing begins here  big arrays! 724 722 iBeg = 0 725 #time0 = time.time()723 time0 = time.time() 726 724 while iBeg < nRef: 727 725 iFin = min(iBeg+blkSize,nRef) … … 862 860 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: 863 861 Flack = 1.2.*parmDict[phfx+'Flack'] 864 #time0 = time.time()862 time0 = time.time() 865 863 #reflection processing begins here  big arrays! 866 864 iBeg = 0 … … 2926 2924 def GetFobsSq(Histograms,Phases,parmDict,calcControls): 2927 2925 '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' 2931 2927 histoList = Histograms.keys() 2932 2928 histoList.sort() 2933 Ka2 = shl = lamRatio = kRatio = None2934 2929 for histogram in histoList: 2935 2930 if 'PWDR' in histogram[:4]: … … 2953 2948 ymb = np.where(ymb,ymb,1.0) 2954 2949 ycmb = np.array(ycyb) 2955 ratio = 1./np.where(ycmb,ycmb/ymb,1.e10) 2950 ratio = 1./np.where(ycmb,ycmb/ymb,1.e10) 2956 2951 refLists = Histogram['Reflection Lists'] 2957 2952 for phase in refLists: … … 2965 2960 phfx = '%d:%d:'%(pId,hId) 2966 2961 refDict = refLists[phase] 2967 sumInt = 0.02968 nExcl = 02969 useMP = G2mp.useMP and len(refDict['RefList']) > 1002970 if useMP: # yes, create a set of initialized Python processes2971 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] *= 12989 nExcl += 12990 elif icod == 2:2991 break2992 if useMP:2993 for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):2994 sumInt += sInt2995 for refl8im,irefl in resList:2996 if refl8im is None:2997 refDict['RefList'][irefl][3+im] *= 12998 nExcl += 12999 else:3000 refDict['RefList'][irefl][8+im] = refl8im3001 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] *= 13016 nExcl += 13017 elif icod == 2:3018 break3019 if useMP:3020 for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):3021 sumInt += sInt3022 for refl8im,irefl in resList:3023 if refl8im is None:3024 refDict['RefList'][irefl][3+im] *= 13025 nExcl += 13026 else:3027 refDict['RefList'][irefl][8+im] = refl8im3028 if parmDict[phfx+'LeBail']:3029 refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]3030 2962 sumFo = 0.0 3031 2963 sumdF = 0.0 3032 2964 sumFosq = 0.0 3033 2965 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 iBegiFin: #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,pos2fmin)) 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 iBegiFin: #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] 3035 3017 Fo = np.sqrt(np.abs(refl[8+im])) 3036 3018 Fc = np.sqrt(np.abs(refl[9]+im)) … … 3051 3033 Histogram = Histograms[histogram] 3052 3034 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 3057 3037 def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): 3058 3038 '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 3063 3042 def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict): 3064 3043 U = parmDict[hfx+'U'] … … 3088 3067 bet = parmDict[hfx+'beta0']+parmDict[hfx+'beta1']/refl[4+im]**4+parmDict[hfx+'betaq']/refl[4+im]**2 3089 3068 return alp,bet 3090 3069 3091 3070 hId = Histogram['hId'] 3092 3071 hfx = ':%d:'%(hId) … … 3096 3075 cw = np.diff(x) 3097 3076 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 3100 3078 if 'C' in calcControls[hfx+'histType']: 3101 3079 shl = max(parmDict[hfx+'SH/L'],0.002) … … 3137 3115 StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict) 3138 3116 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']: 3147 3119 if im: 3148 3120 h,k,l,m = refl[:4] … … 3176 3148 badPeak = True 3177 3149 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 3183 3151 if Ka2: 3184 3152 pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0) # + 360/pi * Dlam/lam * tan(th) … … 3192 3160 elif iBeg > iFin: #bad peak coeff  skip 3193 3161 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']: 3201 3164 h,k,l = refl[:3] 3202 3165 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? 3204 3167 Lorenz = sind(abs(parmDict[hfx+'2theta'])/2)*refl[4+im]**4 #TOF Lorentz correction 3205 3168 # refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict) #apply hydrostatic strain shift 3206 3169 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? 3208 3171 refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) 3209 3172 refl[11+im] *= Vst*Lorenz … … 3228 3191 badPeak = True 3229 3192 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) 3242 3195 if badPeak: 3243 3196 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 3246 3198 return yc,yb 3247 3199 … … 3253 3205 # starttime = time.time() 3254 3206 # 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 3256 3230 # create a list of dependent variables and set up a dictionary to hold their derivatives 3257 3231 dependentVars = G2mv.GetDependentVars() … … 3288 3262 cw = np.append(cw,cw[1]) 3289 3263 Ka2 = False #also for TOF! 3290 dFdvDict = shl = wave = kRatio = lamRatio = None # variables that are passed and must be initialized even when not needed3291 3264 if 'C' in calcControls[hfx+'histType']: 3292 3265 shl = max(parmDict[hfx+'SH/L'],0.002) … … 3298 3271 else: 3299 3272 wave = parmDict[hfx+'Lam'] 3273 #print '#1 getPowderProfileDerv t=',time.time()starttime 3300 3274 for phase in Histogram['Reflection Lists']: 3301 3275 refDict = Histogram['Reflection Lists'][phase] … … 3327 3301 # print 'sfderv time %.3fs'%(time.time()time0) 3328 3302 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 #========================================================================================== 3346 3329 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) 3347 3339 if 'C' in calcControls[hfx+'histType']: #CW powder 3348 3340 Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl) … … 3355 3347 elif not iBegiFin: #peak above high limit  done 3356 3348 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 = iFiniBeg 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,pos2fmin) 3363 iFin2 = np.searchsorted(x,pos2+fmax) 3364 if iBeg2iFin2: 3365 lenBF2 = iFin2iBeg2 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 = iFiniBeg 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+'beta0':[1.0,'bet'],hfx+'beta1':[1./refl[4+im]**4,'bet'], 3415 hfx+'betaq':[1./refl[4+im]**2,'bet'],hfx+'sig0':[1.0,'sig'],hfx+'sig1':[refl[4+im]**2,'sig'], 3416 hfx+'sig2':[refl[4+im]**4,'sig'],hfx+'sigq':[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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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 iFin2iBeg2: 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) 3367 3541 # now process derivatives in constraints 3542 #print '#3 getPowderProfileDerv t=',time.time()starttime 3543 #print timelist,sum(timelist) 3368 3544 dMdv = ma.array(dMdv,mask=np.outer(np.ones(len(varylist)),ma.getmaskarray(x))) #x is a MaskedArray! 3369 3545 G2mv.Dict2Deriv(varylist,depDerivDict,dMdv) … … 3845 4021 M = np.concatenate((M,np.sqrt(pWt)*pVals)) 3846 4022 return M 3847 4023
Note: See TracChangeset
for help on using the changeset viewer.