Changeset 342 for trunk/GSASIIpwd.py
 Timestamp:
 Aug 5, 2011 2:35:43 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIpwd.py
r335 r342 92 92 return plist 93 93 94 def ValuesOut(parmdict, varylist):95 '''Use before call to leastsq to setup list of values for the parameters96 in parmdict, as selected by key in varylist'''97 return [parmdict[key] for key in varylist]98 99 def ValuesIn(parmdict, varylist, values):100 ''' Use after call to leastsq to update the parameter dictionary with101 values corresponding to keys in varylist'''102 parmdict.update(zip(varylist,values))103 104 94 def Transmission(Geometry,Abs,Diam): 105 95 #Calculate sample transmission … … 482 472 fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx') 483 473 484 def getBackground(parmDict,bakType,xdata): 474 def getBackground(pfx,parmDict,bakType,xdata): 475 yb = np.zeros_like(xdata) 485 476 if bakType == 'chebyschev': 486 yb = np.zeros_like(xdata)487 477 iBak = 0 488 478 while True: 489 key = 'Back'+str(iBak)479 key = pfx+'Back:'+str(iBak) 490 480 try: 491 481 yb += parmDict[key]*(xdataxdata[0])**iBak 492 482 iBak += 1 493 483 except KeyError: 494 return yb 495 484 break 485 return yb 486 487 def calcPeakFFT(x,fxns,widths,pos,args): 488 dx = x[1]x[0] 489 z = np.empty([len(fxns),len(x)]) 490 for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)): 491 z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx 492 Z = fft(z) 493 if len(fxns)%2: 494 return ifft(Z.prod(axis=0)).real 495 else: 496 return fftshift(ifft(Z.prod(axis=0))).real 497 498 def getFCJVoigt(pos,intens,sig,gam,shl,xdata): 499 500 DX = xdata[1]xdata[0] 501 widths,fmin,fmax = getWidths(pos,sig,gam,shl) 502 x = np.linspace(posfmin,pos+fmin,1024) 503 if pos > 90: 504 fmin,fmax = [fmax,fmin] 505 dx = x[1]x[0] 506 arg = [pos,shl/10.,dx,] 507 Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0) 508 if len(np.nonzero(Df)[0])>5: 509 fxns = [st.norm,st.cauchy,fcjde] 510 args = [[],[],arg] 511 else: 512 fxns = [st.norm,st.cauchy] 513 args = [[],[]] 514 Df = calcPeakFFT(x,fxns,widths,pos,args) 515 Df /= np.sum(Df) 516 Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0) 517 return intens*Df(xdata)*DX/dx #*10 to get close to old fxn??? 518 496 519 def getPeakProfile(parmDict,xdata,varyList,bakType): 497 520 498 def calcPeakFFT(x,fxns,widths,pos,args): 499 dx = x[1]x[0] 500 z = np.empty([len(fxns),len(x)]) 501 for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)): 502 z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx 503 Z = fft(z) 504 if len(fxns)%2: 505 return ifft(Z.prod(axis=0)).real 506 else: 507 return fftshift(ifft(Z.prod(axis=0))).real 508 509 def getFCJVoigt(pos,intens,sig,gam,shl,xdata): 510 511 DX = xdata[1]xdata[0] 512 widths,fmin,fmax = getWidths(pos,sig,gam,shl) 513 x = np.linspace(posfmin,pos+fmin,1024) 514 if pos > 90: 515 fmin,fmax = [fmax,fmin] 516 dx = x[1]x[0] 517 arg = [pos,shl/10.,dx,] 518 Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0) 519 if len(np.nonzero(Df)[0])>5: 520 fxns = [st.norm,st.cauchy,fcjde] 521 args = [[],[],arg] 522 else: 523 fxns = [st.norm,st.cauchy] 524 args = [[],[]] 525 Df = calcPeakFFT(x,fxns,widths,pos,args) 526 Df /= np.sum(Df) 527 Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0) 528 return intens*Df(xdata)*DX/dx #*10 to get close to old fxn??? 529 530 yb = getBackground(parmDict,bakType,xdata) 521 yb = getBackground('',parmDict,bakType,xdata) 531 522 yc = np.zeros_like(yb) 532 523 U = parmDict['U'] … … 591 582 return widths,fmin,fmax 592 583 593 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data): 584 def Dict2Values(parmdict, varylist): 585 '''Use before call to leastsq to setup list of values for the parameters 586 in parmdict, as selected by key in varylist''' 587 return [parmdict[key] for key in varylist] 588 589 def Values2Dict(parmdict, varylist, values): 590 ''' Use after call to leastsq to update the parameter dictionary with 591 values corresponding to keys in varylist''' 592 parmdict.update(zip(varylist,values)) 593 594 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False): 594 595 595 596 def SetBackgroundParms(Background): 596 597 bakType,bakFlag = Background[:2] 597 598 backVals = Background[3:] 598 backNames = ['Back '+str(i) for i in range(len(backVals))]599 backNames = ['Back:'+str(i) for i in range(len(backVals))] 599 600 if bakFlag: #returns backNames as varyList = backNames 600 601 return bakType,dict(zip(backNames,backVals)),backNames … … 606 607 while True: 607 608 try: 608 bakName = 'Back '+str(iBak)609 bakName = 'Back:'+str(iBak) 609 610 Background[iBak+3] = parmList[bakName] 610 611 iBak += 1 … … 620 621 for i,back in enumerate(Background[3:]): 621 622 ptstr += ptfmt % (back) 622 sigstr += ptfmt % (sigDict['Back '+str(i)])623 sigstr += ptfmt % (sigDict['Back:'+str(i)]) 623 624 print ptstr 624 625 print sigstr … … 675 676 print sigstr 676 677 677 def setPeaksParms(Peaks):678 def SetPeaksParms(Peaks): 678 679 peakNames = [] 679 680 peakVary = [] … … 734 735 return M 735 736 737 Ftol = 0.01 738 if oneCycle: 739 Ftol = 1.0 736 740 x,y,w,yc,yb,yd = data #these are numpy arrays! 737 741 xBeg = np.searchsorted(x,Limits[0]) … … 739 743 bakType,bakDict,bakVary = SetBackgroundParms(Background) 740 744 dataType,insDict,insVary = SetInstParms(Inst) 741 peakDict,peakVary = setPeaksParms(Peaks)745 peakDict,peakVary = SetPeaksParms(Peaks) 742 746 parmDict = {} 743 747 parmDict.update(bakDict) … … 747 751 while True: 748 752 begin = time.time() 749 values = np.array( ValuesOut(parmDict, varyList))753 values = np.array(Dict2Values(parmDict, varyList)) 750 754 if FitPgm == 'LSQ': 751 755 dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0, … … 755 759 dlg.SetPosition(wx.Point(screenSize[2]Size[0]305,screenSize[1]+5)) 756 760 try: 757 result = so.leastsq(errPeakProfile,values, 758 args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg) ,full_output=True)761 result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol, 762 args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg)) 759 763 finally: 760 764 dlg.Destroy() … … 762 766 chisq = np.sum(result[2]['fvec']**2) 763 767 ncyc = int(result[2]['nfev']/len(varyList)) 764 Values In(parmDict, varyList, result[0])768 Values2Dict(parmDict, varyList, result[0]) 765 769 Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100. #to % 766 770 GOF = chisq/(xFinxBeglen(varyList)) … … 786 790 787 791 sigDict = dict(zip(varyList,sig)) 788 yb[xBeg:xFin] = getBackground( parmDict,bakType,x[xBeg:xFin])792 yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin]) 789 793 yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType) 790 794 yd[xBeg:xFin] = y[xBeg:xFin]yc[xBeg:xFin] … … 913 917 msg = 'test ' 914 918 gplot = plotter.add('FCJVoigt, 11BM').gca() 915 gplot.plot(xdata,getBackground( parmDict0,bakType,xdata))919 gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata)) 916 920 gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType)) 917 921 fplot = plotter.add('FCJVoigt, Ka1+2').gca() 918 fplot.plot(xdata,getBackground( parmDict1,bakType,xdata))922 fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata)) 919 923 fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType)) 924 925 def test1(): 926 if NeedTestData: TestData() 920 927 time0 = time.time() 921 928 for i in range(100): 922 929 y = getPeakProfile(parmDict1,xdata,varyList,bakType) 923 print '100+6*Ka12 peaks=3600 peaks',time.time()time0 930 print '100+6*Ka12 peaks=1200 peaks',time.time()time0 931 924 932 925 933
Note: See TracChangeset
for help on using the changeset viewer.