Changeset 353 for trunk/GSASIIpwd.py
 Timestamp:
 Aug 24, 2011 4:07:29 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIpwd.py
r350 r353 492 492 fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx') 493 493 494 def getBackground(pfx,parmDict,bakType,xdata):495 yb = np.zeros_like(xdata)496 if bakType == 'chebyschev':497 iBak = 0498 while True:499 key = pfx+'Back:'+str(iBak)500 try:501 yb += parmDict[key]*(xdataxdata[0])**iBak502 iBak += 1503 except KeyError:504 break505 return yb506 507 494 def getWidths(pos,sig,gam,shl): 508 495 widths = [np.sqrt(sig)/100.,gam/200.] … … 535 522 return intens*Df(xdata)*DX/dx 536 523 524 def getBackground(pfx,parmDict,bakType,xdata): 525 yb = np.zeros_like(xdata) 526 if bakType == 'chebyschev': 527 iBak = 0 528 while True: 529 key = pfx+'Back:'+str(iBak) 530 try: 531 yb += parmDict[key]*(xdataxdata[0])**iBak 532 iBak += 1 533 except KeyError: 534 break 535 return yb 536 537 def getBackgroundDerv(pfx,parmDict,bakType,xdata): 538 dydb = [] 539 if bakType == 'chebyschev': 540 iBak = 0 541 while True: 542 if pfx+'Back:'+str(iBak) in parmDict: 543 dydb.append((xdataxdata[0])**iBak) 544 iBak += 1 545 else: 546 break 547 return dydb 548 537 549 #use old fortran routine 538 def getFCJVoigt3(pos, intens,sig,gam,shl,xdata):550 def getFCJVoigt3(pos,sig,gam,shl,xdata): 539 551 540 552 Df = pyd.pypsvfcj(len(xdata),xdatapos,pos,sig,gam,shl) 541 553 Df /= np.sum(Df) 542 return intens*Df 554 return Df 555 556 def getdFCJVoigt3(pos,sig,gam,shl,xdata): 557 558 Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdatapos,pos,sig,gam,shl) #might have to make these numpy arrays? 559 sumDf = np.sum(Df) 560 return Df/sumDf,dFdp,dFds,dFdg,dFdsh 561 543 562 544 563 def getPeakProfile(parmDict,xdata,varyList,bakType): … … 552 571 X = parmDict['X'] 553 572 Y = parmDict['Y'] 554 shl = max(parmDict['SH/L'],0.00 05)573 shl = max(parmDict['SH/L'],0.002) 555 574 Ka2 = False 556 575 if 'Lam1' in parmDict.keys(): … … 589 608 elif not iBegiFin: #peak above high limit 590 609 return yb+yc 591 yc[iBeg:iFin] += getFCJVoigt3(pos,intens,sig,gam,shl,xdata[iBeg:iFin])610 yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin]) 592 611 if Ka2: 593 612 pos2 = pos+lamRatio*tand(pos/2.0) # + 360/pi * Dlam/lam * tan(th) … … 596 615 iFin = min(lenX,iFin+kdelt) 597 616 if iBegiFin: 598 yc[iBeg:iFin] += getFCJVoigt3(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])617 yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin]) 599 618 iPeak += 1 600 619 except KeyError: #no more peaks to process 601 620 return yb+yc 602 621 622 def getPeakProfileDerv(parmDict,xdata,varyList,bakType): 623 # needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order 624 dMdv = np.zeros(shape=(len(varyList),len(xdata))) 625 if 'Back:0' in varyList: #background derivs are in front if present 626 dMdb = getBackgroundDerv('',parmDict,bakType,xdata) 627 dMdv[0:len(dMdb)] = dMdb 628 629 dx = xdata[1]xdata[0] 630 U = parmDict['U'] 631 V = parmDict['V'] 632 W = parmDict['W'] 633 X = parmDict['X'] 634 Y = parmDict['Y'] 635 shl = max(parmDict['SH/L'],0.002) 636 Ka2 = False 637 if 'Lam1' in parmDict.keys(): 638 Ka2 = True 639 lamRatio = 360*(parmDict['Lam2']parmDict['Lam1'])/(np.pi*parmDict['Lam1']) 640 kRatio = parmDict['I(L2)/I(L1)'] 641 iPeak = 0 642 while True: 643 try: 644 pos = parmDict['pos'+str(iPeak)] 645 intens = parmDict['int'+str(iPeak)] 646 sigName = 'sig'+str(iPeak) 647 tanth = tand(pos/2.0) 648 costh = cosd(pos/2.0) 649 if sigName in varyList: 650 sig = parmDict[sigName] 651 else: 652 sig = U*tanth**2+V*tanth+W 653 dsdU = tanth**2 654 dsdV = tanth 655 dsdW = 1.0 656 sig = max(sig,0.001) #avoid neg sigma 657 gamName = 'gam'+str(iPeak) 658 if gamName in varyList: 659 gam = parmDict[gamName] 660 else: 661 gam = X/costh+Y*tanth 662 dgdX = 1.0/costh 663 dgdY = tanth 664 gam = max(gam,0.001) #avoid neg gamma 665 Wd,fmin,fmax = getWidths(pos,sig,gam,shl) 666 iBeg = np.searchsorted(xdata,posfmin) 667 lenX = len(xdata) 668 if not iBeg: 669 iFin = np.searchsorted(xdata,pos+fmin) 670 elif iBeg == lenX: 671 iFin = iBeg 672 else: 673 iFin = min(lenX,iBeg+int((fmin+fmax)/dx)) 674 if not iBeg+iFin: #peak below low limit 675 iPeak += 1 676 continue 677 elif not iBegiFin: #peak above high limit 678 break 679 dMdpk = np.zeros(shape=(6,len(xdata))) 680 dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin]) 681 for i in range(5): 682 dMdpk[i][iBeg:iFin] = 100.*dx*intens*dMdipk[i] 683 dMdpk[0][iBeg:iFin] = dMdipk[0] 684 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]} 685 if Ka2: 686 pos2 = pos+lamRatio*tand(pos/2.0) # + 360/pi * Dlam/lam * tan(th) 687 kdelt = int((pos2pos)/dx) 688 iBeg = min(lenX,iBeg+kdelt) 689 iFin = min(lenX,iFin+kdelt) 690 if iBegiFin: 691 dMdipk = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin]) 692 for i in range(5): 693 dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk[i] 694 dMdpk[0][iBeg:iFin] += kRatio*dMdipk[0] 695 dMdpk[5][iBeg:iFin] += dMdipk[0] 696 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens} 697 for parmName in ['pos','int','sig','gam']: 698 try: 699 idx = varyList.index(parmName+str(iPeak)) 700 dMdv[idx] = dervDict[parmName] 701 except ValueError: 702 pass 703 if 'U' in varyList: 704 dMdv[varyList.index('U')] += dsdU*dervDict['sig'] 705 if 'V' in varyList: 706 dMdv[varyList.index('V')] += dsdV*dervDict['sig'] 707 if 'W' in varyList: 708 dMdv[varyList.index('W')] += dsdW*dervDict['sig'] 709 if 'X' in varyList: 710 dMdv[varyList.index('X')] += dgdX*dervDict['gam'] 711 if 'Y' in varyList: 712 dMdv[varyList.index('Y')] += dgdY*dervDict['gam'] 713 if 'SH/L' in varyList: 714 dMdv[varyList.index('SH/L')] += dervDict['shl'] #problem here 715 if 'I(L2)/I(L1)' in varyList: 716 dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2'] 717 iPeak += 1 718 except KeyError: #no more peaks to process 719 break 720 return dMdv 721 603 722 def Dict2Values(parmdict, varylist): 604 723 '''Use before call to leastsq to setup list of values for the parameters … … 611 730 parmdict.update(zip(varylist,values)) 612 731 613 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False ):732 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False,controls=None): 614 733 615 734 def SetBackgroundParms(Background): … … 651 770 insVary = [] 652 771 for i,flag in enumerate(insFlags): 653 if flag :772 if flag and insNames[i] in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)']: 654 773 insVary.append(insNames[i]) 655 774 instDict = dict(zip(insNames,insVals)) 656 775 instDict['X'] = max(instDict['X'],0.01) 657 776 instDict['Y'] = max(instDict['Y'],0.01) 658 instDict['SH/L'] = max(instDict['SH/L'],0.00 01)777 instDict['SH/L'] = max(instDict['SH/L'],0.002) 659 778 return dataType,instDict,insVary 660 779 … … 743 862 ptstr += 12*' ' 744 863 print '%s'%(('Peak'+str(i+1)).center(8)),ptstr 864 865 def devPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg): 866 parmdict.update(zip(varylist,values)) 867 return np.sqrt(weights)*getPeakProfileDerv(parmdict,xdata,varylist,bakType) 745 868 746 869 def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg): 747 870 parmdict.update(zip(varylist,values)) 748 M = np.sqrt(weights)*( ydatagetPeakProfile(parmdict,xdata,varylist,bakType))871 M = np.sqrt(weights)*(getPeakProfile(parmdict,xdata,varylist,bakType)ydata) 749 872 Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.) 750 873 if dlg: … … 754 877 return M 755 878 756 Ftol = 0.0001 879 if controls: 880 Ftol = controls['min dM/M'] 881 derivType = controls['deriv type'] 882 else: 883 Ftol = 0.0001 884 derivType = 'analytic' 757 885 if oneCycle: 758 886 Ftol = 1.0 … … 778 906 dlg.SetPosition(wx.Point(screenSize[2]Size[0]305,screenSize[1]+5)) 779 907 try: 780 result = so.leastsq(errPeakProfile,values,full_output=True,epsfcn=1.e8,ftol=Ftol, 781 args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg)) 908 if derivType == 'analytic': 909 result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True, 910 args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg)) 911 ncyc = int(result[2]['nfev']/2) 912 else: 913 result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,epsfcn=1.e8, 914 args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg)) 915 ncyc = int(result[2]['nfev']/len(varyList)) 782 916 finally: 783 917 dlg.Destroy() 784 918 runtime = time.time()begin 785 919 chisq = np.sum(result[2]['fvec']**2) 786 ncyc = int(result[2]['nfev']/len(varyList))787 920 Values2Dict(parmDict, varyList, result[0]) 788 921 Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100. #to % … … 929 1062 'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5, 930 1063 } 1064 global parmDict2 1065 parmDict2 = { 1066 'pos0':5.7,'int0':10.0,'sig0':0.5,'gam0':1.0, 1067 'U':1.,'V':1.,'W':0.1,'X':0.0,'Y':2.,'SH/L':0.004, 1068 'Back0':5.,'Back1':0.02,'Back2':.004, 1069 'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5, 1070 } 931 1071 global varyList 932 1072 varyList = [] … … 949 1089 print '100+6*Ka12 peaks=1200 peaks',time.time()time0 950 1090 1091 def test2(name,delt): 1092 if NeedTestData: TestData() 1093 varyList = [name,] 1094 xdata = np.linspace(5.6,5.8,800) 1095 hplot = plotter.add('derivatives test for '+name).gca() 1096 hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0]) 1097 y0 = getPeakProfile(parmDict2,xdata,varyList,bakType) 1098 parmDict2[name] += delt 1099 y1 = getPeakProfile(parmDict2,xdata,varyList,bakType) 1100 hplot.plot(xdata,(y1y0)/delt,'r+') 1101 951 1102 952 1103 … … 955 1106 global plotter 956 1107 plotter = plot.PlotNotebook() 957 test0() 1108 # test0() 1109 test2('I(L2)/I(L1)',.0001) 958 1110 print "OK" 959 1111 plotter.StartEventLoop()
Note: See TracChangeset
for help on using the changeset viewer.