 Timestamp:
 Mar 5, 2012 12:34:16 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

MPbranch/GSASIIstruct.py
r507 r512 41 41 ''' 42 42 Controls = {'deriv type':'analytic Hessian','max cyc':3, 43 'max Hprocess': 2,44 'max Rprocess': 1,43 'max Hprocess':1, 44 'max Rprocess':2, 45 45 'min dM/M':0.0001,'shift factor':1.} 46 46 file = open(GPXfile,'rb') … … 2000 2000 if pfx+'SHorder' in parmDict: 2001 2001 Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict) 2002 refl[13] = Icorr 2002 # refl[13] = Icorr 2003 return Icorr 2003 2004 2004 2005 def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict): … … 2276 2277 Histogram[phfx+'Nref'] = len(refList) 2277 2278 2278 def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): 2279 2279 def ComputeReflectionProfile(arg): 2280 2280 def GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict): 2281 2281 U = parmDict[hfx+'U'] … … 2290 2290 gam = max(0.001,gam) 2291 2291 return sig,gam 2292 2292 2293 (RefNum,refList,x,Phase,calcControls,wave,G,g,GB,Vst, 2294 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio) = arg 2295 yc = np.zeros_like(x) 2296 iBegO = None 2297 SGData = Phase['General']['SGData'] 2298 shl = max(parmDict[hfx+'SH/L'],0.002) 2299 if 'C' in calcControls[hfx+'histType']: 2300 for refl in refList: 2301 refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position 2302 Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.)) #Lorentz correction 2303 refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift 2304 refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam 2305 refl[13] = Vst*Lorenz * GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13] 2306 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) 2307 iBeg = np.searchsorted(x,refl[5]fmin) 2308 iFin = np.searchsorted(x,refl[5]+fmax) 2309 if not iBeg+iFin: #peak below low limit  skip peak 2310 continue 2311 elif not iBegiFin: #peak above high limit  done 2312 break 2313 yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here 2314 if iBegO is None: 2315 iBegO = iBeg 2316 iFinO = iFin 2317 else: 2318 iBegO = min(iBegO,iBeg) 2319 iFinO = max(iFinO,iFin) 2320 if Ka2: 2321 pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) 2322 Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl) 2323 iBeg2 = np.searchsorted(x,pos2fmin) 2324 iFin2 = np.searchsorted(x,pos2+fmax) 2325 if not iBeg2+iFin2: #peak below low limit  skip peak 2326 continue 2327 elif not iBeg2iFin2: #peak above high limit  done 2328 continue 2329 yc[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here 2330 iBegO = min(iBegO,iBeg2) 2331 iFinO = max(iFinO,iFin2) 2332 elif 'T' in calcControls[hfx+'histType']: 2333 print 'TOF Undefined at present' 2334 raise Exception #no TOF yet 2335 2336 if iBegO is None: return None # no reflections were computed 2337 return RefNum,refList,iBegO,iFinO,yc[iBegO:iFinO] 2338 2339 def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): 2340 mpPool = mp.Pool() 2293 2341 hId = Histogram['hId'] 2294 2342 hfx = ':%d:'%(hId) … … 2299 2347 if 'C' in calcControls[hfx+'histType']: 2300 2348 shl = max(parmDict[hfx+'SH/L'],0.002) 2301 Ka2 = False2302 2349 if hfx+'Lam1' in parmDict.keys(): 2303 2350 wave = parmDict[hfx+'Lam1'] … … 2307 2354 else: 2308 2355 wave = parmDict[hfx+'Lam'] 2356 Ka2 = False 2357 lamRatio = None 2358 kRatio = None 2309 2359 else: 2310 2360 print 'TOF Undefined at present' … … 2324 2374 if 'Pawley' not in Phase['General']['Type']: 2325 2375 refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict) 2326 for refl in refList:2327 if 'C' in calcControls[hfx+'histType']:2376 else: 2377 for refl in refList: 2328 2378 h,k,l = refl[:3] 2329 refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position 2330 Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.)) #Lorentz correction 2331 refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift 2332 refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam 2333 GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13] 2334 refl[13] *= Vst*Lorenz 2335 if 'Pawley' in Phase['General']['Type']: 2336 try: 2337 refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]) 2338 except KeyError: 2339 # print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) 2340 continue 2341 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) 2342 iBeg = np.searchsorted(x,refl[5]fmin) 2343 iFin = np.searchsorted(x,refl[5]+fmax) 2344 if not iBeg+iFin: #peak below low limit  skip peak 2379 try: 2380 refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]) 2381 except KeyError: 2382 #print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) 2345 2383 continue 2346 elif not iBegiFin: #peak above high limit  done 2347 break 2348 yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here 2349 if Ka2: 2350 pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) 2351 Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl) 2352 iBeg = np.searchsorted(x,pos2fmin) 2353 iFin = np.searchsorted(x,pos2+fmax) 2354 if not iBeg+iFin: #peak below low limit  skip peak 2355 continue 2356 elif not iBegiFin: #peak above high limit  done 2357 return yc,yb 2358 yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here 2359 elif 'T' in calcControls[hfx+'histType']: 2360 print 'TOF Undefined at present' 2361 raise Exception #no TOF yet 2384 2385 chunk = 200 2386 if (calcControls['max Hprocess'] > 1 or 2387 calcControls['max Rprocess'] == 1 or 2388 mp.cpu_count() < 2 or len(refList) < 2 * chunk): 2389 print 'single thread function',hId,len(refList) 2390 res = ComputeReflectionProfile( 2391 (0,refList,x,Phase,calcControls,wave,G,g,GB,Vst, 2392 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio)) 2393 if not res: continue 2394 RefNum,refls,iBeg,iFin,yct = res 2395 yc[iBeg:iFin] += yct 2396 else: 2397 print 'multiprocess function',hId,len(refList) 2398 argList = [] 2399 for i in range(0,len(refList),chunk): 2400 argList.append( 2401 (i,refList[i:i+chunk],x,Phase,calcControls,wave,G,g,GB,Vst, 2402 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio)) 2403 for res in mpPool.imap_unordered(ComputeReflectionProfile,argList): 2404 if not res: continue 2405 RefNum,refls,iBeg,iFin,yct = res 2406 yc[iBeg:iFin] += yct 2407 for i,refl in enumerate(refls): 2408 refList[RefNum+i] = refl 2409 mpPool.close() 2410 mpPool.join() 2362 2411 return yc,yb 2363 2364 def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): 2365 2412 2413 def ComputeReflectionDerivative(arg): 2366 2414 def cellVaryDerv(pfx,SGData,dpdA): 2367 2415 if SGData['SGLaue'] in ['1',]: … … 2388 2436 # return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]] 2389 2437 return [[pfx+'A0',dpdA[0]]] 2438 2439 (RefNum,refList,x,dependentVars,varylist,dFdvDict, 2440 Phase,calcControls,wave,G,g,GB,A, 2441 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio) = arg 2442 dMdv = np.zeros(shape=(len(varylist),len(x))) 2443 depDerivDict = {} 2444 for j in dependentVars: 2445 depDerivDict[j] = np.zeros(shape=(len(x))) 2446 SGData = Phase['General']['SGData'] 2447 shl = max(parmDict[hfx+'SH/L'],0.002) 2448 dx = x[1]x[0] 2449 iBegO = None 2450 for irefS,refl in enumerate(refList): 2451 iref = irefS + RefNum 2452 if 'C' in calcControls[hfx+'histType']: #CW powder 2453 dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) 2454 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) 2455 iBeg = np.searchsorted(x,refl[5]fmin) 2456 iFin = np.searchsorted(x,refl[5]+fmax) 2457 if iBegO is None: 2458 iBegO = iBeg 2459 iFinO = iFin 2460 else: 2461 iBegO = min(iBegO,iBeg) 2462 iFinO = max(iFinO,iFin) 2463 if not iBeg+iFin: #peak below low limit  skip peak 2464 continue 2465 #return None 2466 elif not iBegiFin: #peak above high limit  done 2467 break 2468 #return None 2469 pos = refl[5] 2470 tanth = tand(pos/2.0) 2471 costh = cosd(pos/2.0) 2472 lenBF = iFiniBeg 2473 dMdpk = np.zeros(shape=(6,lenBF)) 2474 dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) 2475 for i in range(1,5): 2476 dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i] 2477 dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0] 2478 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])} 2479 if Ka2: 2480 pos2 = refl[5]+lamRatio*tanth # + 360/pi * Dlam/lam * tan(th) 2481 kdelt = int((pos2refl[5])/dx) 2482 iBeg2 = min(lenX,iBeg+kdelt) 2483 iFin2 = min(lenX,iFin+kdelt) 2484 if iBeg2iFin2: 2485 iBegO = min(iBegO,iBeg2) 2486 iFinO = max(iFinO,iFin2) 2487 lenBF2 = iFin2iBeg2 2488 dMdpk2 = np.zeros(shape=(6,lenBF2)) 2489 dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) 2490 for i in range(1,5): 2491 dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i] 2492 dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0] 2493 dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0] 2494 dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]} 2495 if 'Pawley' in Phase['General']['Type']: 2496 try: 2497 idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%refl[:3]])) 2498 dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9] 2499 if Ka2: 2500 dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9] 2501 # Assuming Pawley variables not in constraints 2502 except ValueError: 2503 pass 2504 dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict) 2505 names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'], 2506 hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'], 2507 hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'], 2508 hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'], 2509 hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'], 2510 hfx+'DisplaceY':[dpdY,'pos'],} 2511 for name in names: 2512 item = names[name] 2513 if name in varylist: 2514 dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]] 2515 if Ka2: 2516 dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]] 2517 elif name in dependentVars: 2518 if Ka2: 2519 depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]] 2520 depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]] 2521 for iPO in dIdPO: 2522 if iPO in varylist: 2523 dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] 2524 if Ka2: 2525 dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int'] 2526 elif iPO in dependentVars: 2527 depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] 2528 if Ka2: 2529 depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int'] 2530 for i,name in enumerate(['omega','chi','phi']): 2531 aname = pfx+'SH '+name 2532 if aname in varylist: 2533 dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int'] 2534 if Ka2: 2535 dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int'] 2536 elif aname in dependentVars: 2537 depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int'] 2538 if Ka2: 2539 depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int'] 2540 for iSH in dFdODF: 2541 if iSH in varylist: 2542 dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] 2543 if Ka2: 2544 dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int'] 2545 elif iSH in dependentVars: 2546 depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] 2547 if Ka2: 2548 depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int'] 2549 cellDervNames = cellVaryDerv(pfx,SGData,dpdA) 2550 for name,dpdA in cellDervNames: 2551 if name in varylist: 2552 dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos'] 2553 if Ka2: 2554 dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos'] 2555 elif name in dependentVars: 2556 depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos'] 2557 if Ka2: 2558 depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos'] 2559 dDijDict = GetHStrainShiftDerv(refl,SGData,phfx) 2560 for name in dDijDict: 2561 if name in varylist: 2562 dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos'] 2563 if Ka2: 2564 dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos'] 2565 elif name in dependentVars: 2566 depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos'] 2567 if Ka2: 2568 depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos'] 2569 gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict) 2570 for name in gamDict: 2571 if name in varylist: 2572 dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam'] 2573 if Ka2: 2574 dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam'] 2575 elif name in dependentVars: 2576 depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam'] 2577 if Ka2: 2578 depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam'] 2579 2580 elif 'T' in calcControls[hfx+'histType']: 2581 print 'TOF Undefined at present' 2582 raise Exception #no TOF yet 2583 #do atom derivatives  for F,X & U so far 2584 corr = dervDict['int']/refl[9] 2585 if Ka2: 2586 corr2 = dervDict2['int']/refl[9] 2587 for name in varylist+dependentVars: 2588 try: 2589 aname = name.split(pfx)[1][:2] 2590 if aname not in ['Af','dA','AU']: continue # skip anything not an atom param 2591 except IndexError: 2592 continue 2593 if name in varylist: 2594 dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr 2595 if Ka2: 2596 dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2 2597 elif name in dependentVars: 2598 depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr 2599 if Ka2: 2600 depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2 2601 # end loop over reflections 2602 return depDerivDict,dMdv[:,iBegO:iFinO],iBegO,iFinO 2603 2604 def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): 2390 2605 # create a list of dependent variables and set up a dictionary to hold their derivatives 2391 2606 dependentVars = G2mv.GetDependentVars() … … 2413 2628 dx = x[1]x[0] 2414 2629 shl = max(parmDict[hfx+'SH/L'],0.002) 2415 Ka2 = False2416 2630 if hfx+'Lam1' in parmDict.keys(): 2417 2631 wave = parmDict[hfx+'Lam1'] … … 2421 2635 else: 2422 2636 wave = parmDict[hfx+'Lam'] 2637 Ka2 = False 2638 lamRatio = None 2639 kRatio = None 2423 2640 else: 2424 2641 print 'TOF Undefined at present' … … 2436 2653 if 'Pawley' not in Phase['General']['Type']: 2437 2654 dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict) 2438 for iref,refl in enumerate(refList): 2439 if 'C' in calcControls[hfx+'histType']: #CW powder 2440 h,k,l = refl[:3] 2441 dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) 2442 if 'Pawley' in Phase['General']['Type']: 2443 try: 2444 refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]) 2445 except KeyError: 2446 # print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) 2447 continue 2448 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) 2449 iBeg = np.searchsorted(x,refl[5]fmin) 2450 iFin = np.searchsorted(x,refl[5]+fmax) 2451 if not iBeg+iFin: #peak below low limit  skip peak 2655 else: 2656 for iref,refl in enumerate(refList): 2657 try: 2658 refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%refl[:3]])]) 2659 except KeyError: 2660 #print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%refl[:3] 2452 2661 continue 2453 elif not iBegiFin: #peak above high limit  done 2454 break 2455 pos = refl[5] 2456 tanth = tand(pos/2.0) 2457 costh = cosd(pos/2.0) 2458 lenBF = iFiniBeg 2459 dMdpk = np.zeros(shape=(6,lenBF)) 2460 dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) 2461 for i in range(1,5): 2462 dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i] 2463 dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0] 2464 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])} 2465 if Ka2: 2466 pos2 = refl[5]+lamRatio*tanth # + 360/pi * Dlam/lam * tan(th) 2467 kdelt = int((pos2refl[5])/dx) 2468 iBeg2 = min(lenX,iBeg+kdelt) 2469 iFin2 = min(lenX,iFin+kdelt) 2470 if iBeg2iFin2: 2471 lenBF2 = iFin2iBeg2 2472 dMdpk2 = np.zeros(shape=(6,lenBF2)) 2473 dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) 2474 for i in range(1,5): 2475 dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i] 2476 dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0] 2477 dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0] 2478 dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]} 2479 if 'Pawley' in Phase['General']['Type']: 2480 try: 2481 idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])) 2482 dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9] 2483 if Ka2: 2484 dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9] 2485 # Assuming Pawley variables not in constraints 2486 except ValueError: 2487 pass 2488 dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict) 2489 names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'], 2490 hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'], 2491 hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'], 2492 hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'], 2493 hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'], 2494 hfx+'DisplaceY':[dpdY,'pos'],} 2495 for name in names: 2496 item = names[name] 2497 if name in varylist: 2498 dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]] 2499 if Ka2: 2500 dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]] 2501 elif name in dependentVars: 2502 if Ka2: 2503 depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]] 2504 depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]] 2505 for iPO in dIdPO: 2506 if iPO in varylist: 2507 dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] 2508 if Ka2: 2509 dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int'] 2510 elif iPO in dependentVars: 2511 depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] 2512 if Ka2: 2513 depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int'] 2514 for i,name in enumerate(['omega','chi','phi']): 2515 aname = pfx+'SH '+name 2516 if aname in varylist: 2517 dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int'] 2518 if Ka2: 2519 dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int'] 2520 elif aname in dependentVars: 2521 depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int'] 2522 if Ka2: 2523 depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int'] 2524 for iSH in dFdODF: 2525 if iSH in varylist: 2526 dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] 2527 if Ka2: 2528 dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int'] 2529 elif iSH in dependentVars: 2530 depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] 2531 if Ka2: 2532 depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int'] 2533 cellDervNames = cellVaryDerv(pfx,SGData,dpdA) 2534 for name,dpdA in cellDervNames: 2535 if name in varylist: 2536 dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos'] 2537 if Ka2: 2538 dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos'] 2539 elif name in dependentVars: 2540 depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos'] 2541 if Ka2: 2542 depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos'] 2543 dDijDict = GetHStrainShiftDerv(refl,SGData,phfx) 2544 for name in dDijDict: 2545 if name in varylist: 2546 dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos'] 2547 if Ka2: 2548 dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos'] 2549 elif name in dependentVars: 2550 depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos'] 2551 if Ka2: 2552 depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos'] 2553 gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict) 2554 for name in gamDict: 2555 if name in varylist: 2556 dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam'] 2557 if Ka2: 2558 dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam'] 2559 elif name in dependentVars: 2560 depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam'] 2561 if Ka2: 2562 depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam'] 2563 2564 elif 'T' in calcControls[hfx+'histType']: 2565 print 'TOF Undefined at present' 2566 raise Exception #no TOF yet 2567 #do atom derivatives  for F,X & U so far 2568 corr = dervDict['int']/refl[9] 2569 if Ka2: 2570 corr2 = dervDict2['int']/refl[9] 2571 for name in varylist+dependentVars: 2572 try: 2573 aname = name.split(pfx)[1][:2] 2574 if aname not in ['Af','dA','AU']: continue # skip anything not an atom param 2575 except IndexError: 2576 continue 2577 if name in varylist: 2578 dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr 2579 if Ka2: 2580 dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2 2581 elif name in dependentVars: 2582 depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr 2583 if Ka2: 2584 depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2 2662 chunk = 200 2663 if (calcControls['max Hprocess'] > 1 or 2664 calcControls['max Rprocess'] == 1 or 2665 mp.cpu_count() < 2 or len(refList) < 2 * chunk): 2666 print 'single thread derivs' 2667 TdepDerivDict, TdMdv, iBegO, iFinO = ComputeReflectionDerivative( 2668 (0,refList,x,dependentVars,varylist,dFdvDict, 2669 Phase,calcControls,wave,G,g,GB,A, 2670 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio)) 2671 dMdv[:,iBegO:iFinO] += TdMdv 2672 for j in dependentVars: 2673 depDerivDict[j] += TdepDerivDict[j] 2674 else: 2675 print 'multiprocess derivs' 2676 argList = [] 2677 for iref in range(0,len(refList),chunk): 2678 argList.append( 2679 (iref,refList[iref:iref+chunk],x,dependentVars,varylist,dFdvDict, 2680 Phase,calcControls,wave,G,g,GB,A, 2681 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio)) 2682 mpPool = mp.Pool() 2683 for res in mpPool.imap_unordered(ComputeReflectionDerivative,argList): 2684 if res is None: continue 2685 TdepDerivDict, TdMdv, iBegO, iFinO = res 2686 dMdv[:,iBegO:iFinO] += TdMdv 2687 for j in dependentVars: 2688 depDerivDict[j] += TdepDerivDict[j] 2689 2585 2690 # now process derivatives in constraints 2586 2691 G2mv.Dict2Deriv(varylist,depDerivDict,dMdv) … … 2675 2780 Histogram,parmdict,varylist,Phases, 2676 2781 calcControls,pawleyLookup]) 2677 if MaxProcess > 1: 2678 mpPool = mp.Pool(processes=MaxProcess) 2679 #results = mpPool.map(ComputePowderHessian,argList) 2680 #for Vec,Hess in results: 2782 if MaxProcess > 1: 2783 mpPool = mp.Pool() 2681 2784 for Vec,Hess in mpPool.imap_unordered(ComputePowderHessian,argList): 2682 2785 VecSum += Vec 2683 2786 HessSum += Hess 2787 mpPool.close() 2788 mpPool.join() 2684 2789 else: 2685 2790 for arg in argList: … … 2767 2872 ) 2768 2873 if MaxProcess > 1: 2769 mpPool = mp.Pool( processes=MaxProcess)2874 mpPool = mp.Pool() 2770 2875 for (xB,xF,ycSect,ybSect,RL,histkey 2771 2876 ) in mpPool.imap_unordered(ComputePowderProfile,argList): … … 2789 2894 Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.) 2790 2895 M = np.concatenate((M,wdy)) 2896 mpPool.close() 2897 mpPool.join() 2791 2898 else: 2792 2899 for arg in argList: … … 3313 3420 print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]) 3314 3421 3315 3316 3422 def main(): 3317 3423 arg = sys.argv
Note: See TracChangeset
for help on using the changeset viewer.