Changeset 296 for trunk/GSASIIpwd.py
 Timestamp:
 Jun 9, 2011 4:02:38 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIpwd.py
r292 r296 1 #/usr/bin/env python 2 # * coding: utf8 * 1 3 #GSASII powder calculation module 2 4 ########### SVN repository information ################### … … 11 13 import wx 12 14 import time 15 13 16 import numpy as np 14 17 import scipy as sp 15 18 import numpy.linalg as nl 19 from numpy.fft import ifft, fft, fftshift 16 20 import scipy.interpolate as si 17 21 import scipy.stats as st 22 import scipy.optimize as so 23 18 24 import GSASIIpath 19 import pypowder as pyp #assumes path has been amended to include correctr bin directory 25 #import pypowder as pyp 20 26 import GSASIIplot as G2plt 21 27 import GSASIIlattice as G2lat 22 28 import GSASIIElem as G2elem 23 29 import GSASIIgrid as G2gd 30 import GSASIIIO as G2IO 24 31 25 32 # trig functions in degrees … … 43 50 npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave) 44 51 45 np.seterr(divide='ignore') #this is for the FCJ functions 46 47 #Peak shape definitions 48 # FingerCox_Jephcoat D(2phi,2th) function for S/L = H/L 49 50 class fcjde_gen(st.rv_continuous): 51 52 """ 53 FingerCoxJephcoat D(2phi,2th) function for S/L = H/L 54 Ref: J. Appl. Cryst. (1994) 27, 892900. 55 Parameters 56  57 x: array like 2theta steps (1 to 1) 58 t: 2theta position of peak 59 s: sum(S/L,H/L); S: sample height, H: detector opening, 60 L: sample to detector opening distance 61 dx: 2theta step size in deg 62 Result for fcj.pdf 63  64 if t < 90: 65 X = dx*x+t 66 else: 67 X = 180.dx*xt 68 if X < t & s = S/L+H/L: 69 fcj.pdf = [1/sqrt({cos(x)**2/cos(t)**2}1)  1/s]/cos(X) 70 if X >= t: 71 fcj.pdf = 0 72 """ 73 def _pdf(self,x,t,s,dx): 74 T = np.where(t<=90.,dx*x+t,180.dx*xt) 75 ax = npcosd(T)**2 76 bx = npcosd(t)**2 77 bx = np.where(ax>bx,bx,ax) 78 fx = np.where(ax>bx,(np.sqrt(bx/(axbx))1./s)/np.sqrt(ax),0.0) 79 fx = np.where(fx > 0,fx,0.0) 80 return fx 81 # def _cdf(self, x): 82 # def _ppf(self, q): 83 # def _sf(self, x): 84 # def _isf(self, q): 85 # def _stats(self): 86 # def _entropy(self): 87 88 fcjde = fcjde_gen(name='fcjde') 89 90 91 # FingerCox_Jephcoat D(2phi,2th) function for S/L != H/L 92 93 class fcjd_gen(st.rv_continuous): 94 """ 95 FingerCoxJephcoat D(2phi,2th) function for S/L != H/L 96 Ref: J. Appl. Cryst. (1994) 27, 892900. 97 Parameters 98  99 x: array like 2theta step numbers (1 to 1) 100 t: 2theta position of peak 101 h: detector opening height/sample to detector opening distance 102 s: sample height/sample to detector opening distance 103 dx: 2theta step size in deg 104 Result for fcj2.pdf 105  106 if t < 90: 107 X = dx*x+t 108 else: 109 X = 180.dx*xt 110 infl = acos(cos(t)*sqrt((hs)**2+1)) 111 if X < infl: 112 fcj.pdf = [1/sqrt({cos(X)**2/cos(t)**2}1)  1/shl]/cos(X) 113 for X < t & s = S/L+H/L 114 fcj.pdf(x,t,s) = 0 115 for 2phi >= 2th 116 """ 117 def _pdf(self,x,t,h,s,dx): 118 T = np.where(t<=90.,dx*x+t,180.dx*xt) 119 a = np.abs(npcosd(t))*(np.sqrt((hs)**2+1.)) 120 infl = np.where((a <= 1.),np.abs(npacosd(a)),T) 121 ax = npcosd(T)**2 122 bx = npcosd(t)**2 123 bx = np.where(ax>bx,bx,ax) 124 H = np.where(ax>bx,np.sqrt((ax/bx)1.),0.0) 125 W1 = h+sH 126 W2 = np.where ((h > s),2.*s,2.*h) 127 fx = 2.*h*np.sqrt((ax/bx)1.)*np.sqrt(ax) 128 fx = np.where(fx>0.0,1./fx,0.0) 129 fx = np.where((T < infl),fx*W1,fx*W2) 130 return np.where((fx > 0.),fx,0.0) 131 # def _cdf(self, x): 132 # def _ppf(self, q): 133 # def _sf(self, x): 134 # def _isf(self, q): 135 # def _stats(self): 136 # def _entropy(self): 137 138 fcjd = fcjd_gen(name='fcjd') 139 140 # FingerCox_Jephcoat D(2phi,2th) function for S/L != H/L using sum & difference 141 142 class fcjdsd_gen(st.rv_continuous): 143 """ 144 FingerCoxJephcoat D(2phi,2th) function for S/L != H/L using sum & difference 145 Ref: J. Appl. Cryst. (1994) 27, 892900. 146 Parameters 147  148 x: array like 2theta positions 149 t: 2theta position of peak 150 s: sum(S/L,H/L); S: sample height, H: detector opening, 151 L: sample to detector opening distance 152 d: difference(S/L,H/L) 153 dx: 2theta step size in deg 154 Result for fcj2.pdf 155  156 if t < 90: 157 X = dx*x+t 158 else: 159 X = 180.dx*xt 160 161 fcj.pdf(x,t,s,d) = [1/sqrt({cos(X)**2/cos(t)**2}1)  1/shl]/cos(2phi) 162 163 for 2phi < 2tth & shl = S/L+H/L 164 165 fcj.pdf(x,tth,shl) = 0 166 167 for 2phi >= 2th 168 """ 169 def _argcheck(self,t,s,d,dx): 170 return (t > 0)&(s > 0)&(abs(d) < s) 171 def _pdf(self,x,t,s,d,dx): 172 T = np.where(t<=90.,dx*x+t,180.dx*xt) 173 a = np.abs(npcosd(t))*np.sqrt(d**2+1.) 174 infl = np.where((a < 1.),np.abs(npacosd(a)),T) 175 ax = npcosd(T)**2 176 bx = npcosd(t)**2 177 bx = np.where(ax>bx,bx,ax) 178 H = np.where(ax>bx,np.sqrt((ax/bx)1.),0.0) 179 W1 = sH 180 W2 = np.where ((d > 0),sd,s+d) 181 fx = np.where(ax>bx,1./((s+d)*np.sqrt((ax/bx)1.)*np.sqrt(ax)),0.0) 182 fx = np.where((T < infl),fx*W1,fx*W2) 183 return np.where((fx > 0.),fx,0.0) 184 # def _cdf(self, x): 185 # def _ppf(self, q): 186 # def _sf(self, x): 187 # def _isf(self, q): 188 # def _stats(self): 189 # def _entropy(self): 190 191 fcjdsd = fcjdsd_gen(name='fcjdsd') 192 52 193 53 def factorize(num): 194 54 ''' Provide prime number factors for integer num … … 233 93 return plist 234 94 95 def ValuesOut(parmdict, varylist): 96 '''Use before call to leastsq to setup list of values for the parameters 97 in parmdict, as selected by key in varylist''' 98 return [parmdict[key] for key in varylist] 99 100 def ValuesIn(parmdict, varylist, values): 101 ''' Use after call to leastsq to update the parameter dictionary with 102 values corresponding to keys in varylist''' 103 parmdict.update(zip(varylist,values)) 104 235 105 def Transmission(Geometry,Abs,Diam): 236 106 #Calculate sample transmission … … 547 417 Amu += el['FormulaNo']*el['mu'] 548 418 549 550 def ValEsd(value,esd=0,nTZ=False): #NOT complete  don't use 551 # returns value(esd) string; nTZ=True for no trailing zeros 552 # use esd < 0 for level of precision shown e.g. esd=0.01 gives 2 places beyond decimal 553 #get the 2 significant digits in the esd 554 edig = lambda esd: int(round(10**(math.log10(esd) % 1+1))) 555 #get the number of digits to represent them 556 epl = lambda esd: 2+int(1.545math.log10(10*edig(esd))) 557 558 mdec = lambda esd: int(math.log10(abs(esd))) 559 ndec = lambda esd: int(1.545math.log10(abs(esd))) 560 if esd > 0: 561 fmt = '"%.'+str(ndec(esd))+'f(%d)"' 562 print fmt,ndec(esd),esd*10**(mdec(esd)+1) 563 return fmt%(value,int(esd*10**(mdec(esd)+1))) 564 elif esd < 0: 565 return str(round(value,mdec(esd))) 419 #GSASII peak fitting routine: Finger, Cox & Jephcoat model 420 421 np.seterr(divide='ignore') 422 423 # Voigt function = convolution (norm,cauchy) 424 425 class voigt_gen(st.rv_continuous): 426 ''' 427 Voigt function 428 Parameters 429  430 ''' 431 def _pdf(self,x,s,g): 432 z = np.empty([2,len(x)]) 433 z[0] = st.norm.pdf(x,scale=s) 434 z[1] = st.cauchy.pdf(x,scale=g) 435 Z = fft(z) 436 return fftshift(ifft(Z.prod(axis=0))).real 437 # def _cdf(self, x): 438 # def _ppf(self, q): 439 # def _sf(self, x): 440 # def _isf(self, q): 441 # def _stats(self): 442 # def _entropy(self): 443 444 voigt = voigt_gen(name='voigt',shapes='ts,g') 445 446 # FingerCox_Jephcoat D(2phi,2th) function for S/L = H/L 447 448 class fcjde_gen(st.rv_continuous): 449 """ 450 FingerCoxJephcoat D(2phi,2th) function for S/L = H/L 451 Ref: J. Appl. Cryst. (1994) 27, 892900. 452 Parameters 453  454 x: array 1 to 1 455 t: 2theta position of peak 456 s: sum(S/L,H/L); S: sample height, H: detector opening, 457 L: sample to detector opening distance 458 dx: 2theta step size in deg 459 Result for fcj.pdf 460  461 if x < t & s = S/L+H/L: 462 fcj.pdf = [1/sqrt({cos(x)**2/cos(t)**2}1)  1/s]/cos(x) 463 if x >= t: 464 fcj.pdf = 0 465 """ 466 def _pdf(self,x,t,s,dx): 467 # T = np.where(t<=90.,dx*x+t,180.dx*xt) 468 T = dx*x+t 469 ax = npcosd(T)**2 470 bx = npcosd(t)**2 471 bx = np.where(ax>bx,bx,ax) 472 fx = np.where(ax>bx,(np.sqrt(bx/(axbx))1./s)/np.sqrt(ax),0.0) 473 fx = np.where(fx > 0.,fx,0.0) 474 return fx 475 # def _cdf(self, x): 476 # def _ppf(self, q): 477 # def _sf(self, x): 478 # def _isf(self, q): 479 # def _stats(self): 480 # def _entropy(self): 481 482 fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx') 483 484 def getBackground(parmDict,bakType,xdata): 485 if bakType == 'chebyschev': 486 yb = np.zeros_like(xdata) 487 iBak = 0 488 while True: 489 key = 'Back'+str(iBak) 490 try: 491 yb += parmDict[key]*(xdataxdata[0])**iBak 492 iBak += 1 493 except KeyError: 494 return yb 495 496 def getPeakProfile(parmDict,xdata,varyList,bakType): 497 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) 531 yc = np.zeros_like(yb) 532 U = parmDict['U'] 533 V = parmDict['V'] 534 W = parmDict['W'] 535 X = parmDict['X'] 536 Y = parmDict['Y'] 537 shl = max(parmDict['SH/L'],0.001) 538 Ka2 = False 539 if 'Lam1' in parmDict.keys(): 540 Ka2 = True 541 lamRatio = parmDict['Lam2']/parmDict['Lam1'] 542 kRatio = parmDict['I(L2)/I(L1)'] 543 iPeak = 0 544 while True: 545 try: 546 pos = parmDict['pos'+str(iPeak)] 547 intens = parmDict['int'+str(iPeak)] 548 sigName = 'sig'+str(iPeak) 549 if sigName in varyList: 550 sig = parmDict[sigName] 551 else: 552 sig = U*tand(pos/2.0)**2+V*tand(pos/2.0)+W 553 sig = max(sig,0.001) #avoid neg sigma 554 gamName = 'gam'+str(iPeak) 555 if gamName in varyList: 556 gam = parmDict[gamName] 557 else: 558 gam = X/cosd(pos/2.0)+Y*tand(pos/2.0) 559 gam = max(gam,0.01) #avoid neg gamma 560 Wd,fmin,fmax = getWidths(pos,sig,gam,shl) 561 if pos > 90: 562 fmin,fmax = [fmax,fmin] 563 iBeg = np.searchsorted(xdata,posfmin) 564 iFin = np.searchsorted(xdata,pos+fmax) 565 if not iBeg+iFin: #peak below low limit 566 iPeak += 1 567 continue 568 elif not iBegiFin: #peak above high limit 569 return yb+yc 570 yc[iBeg:iFin] += getFCJVoigt(pos,intens,sig,gam,shl,xdata[iBeg:iFin]) 571 if Ka2: 572 pos2 = 2.0*asind(lamRatio*sind(pos/2.0)) 573 Wd,fmin,fmax = getWidths(pos2,sig,gam,shl) 574 if pos > 90: 575 fmin,fmax = [fmax,fmin] 576 iBeg = np.searchsorted(xdata,pos2fmin) 577 iFin = np.searchsorted(xdata,pos2+fmax) 578 yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin]) 579 iPeak += 1 580 except KeyError: #no more peaks to process 581 return yb+yc 582 583 def getWidths(pos,sig,gam,shl): 584 if shl: 585 widths = [np.sqrt(sig)/100.,gam/200.,.1] 566 586 else: 567 text = "%F"%(value) 568 if nTZ: 569 return text.rstrip('0') 570 else: 571 return text 572 573 574 #GSASII peak fitting routine: Thompson, Cox & Hastings; Finger, Cox & Jephcoat model 575 576 def DoPeakFit(peaks,background,limits,inst,data): 577 578 def backgroundPrint(background,sigback): 579 if background[1]: 580 print 'Background coefficients for',background[0],'function' 587 widths = [np.sqrt(sig)/100.,gam/200.] 588 fwhm = 2.355*widths[0]+2.*widths[1] 589 fmin = 10.*(fwhm+shl*abs(npcosd(pos))) 590 fmax = 15.0*fwhm 591 return widths,fmin,fmax 592 593 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data): 594 595 def SetBackgroundParms(Background): 596 bakType,bakFlag = Background[:2] 597 backVals = Background[3:] 598 backNames = ['Back'+str(i) for i in range(len(backVals))] 599 if bakFlag: #returns backNames as varyList = backNames 600 return bakType,dict(zip(backNames,backVals)),backNames 601 else: #no background varied; varyList = [] 602 return bakType,dict(zip(backNames,backVals)),[] 603 604 def GetBackgroundParms(parmList,Background): 605 iBak = 0 606 while True: 607 try: 608 bakName = 'Back'+str(iBak) 609 Background[iBak+3] = parmList[bakName] 610 iBak += 1 611 except KeyError: 612 break 613 614 def BackgroundPrint(Background,sigDict): 615 if Background[1]: 616 print 'Background coefficients for',Background[0],'function' 581 617 ptfmt = "%12.5f" 582 618 ptstr = 'values:' 583 619 sigstr = 'esds :' 584 for i,back in enumerate( background[3:]):620 for i,back in enumerate(Background[3:]): 585 621 ptstr += ptfmt % (back) 586 sigstr += ptfmt % (sig back[i+3])622 sigstr += ptfmt % (sigDict['Back'+str(i)]) 587 623 print ptstr 588 624 print sigstr 589 625 else: 590 626 print 'Background not refined' 591 592 def instPrint(instVal,siginst,insLabels): 627 628 def SetInstParms(Inst): 629 insVals,insFlags,insNames = Inst[1:4] 630 dataType = insVals[0] 631 insVary = [] 632 for i,flag in enumerate(insFlags): 633 if flag: 634 insVary.append(insNames[i]) 635 instDict = dict(zip(insNames,insVals)) 636 instDict['X'] = max(instDict['X'],0.1) 637 instDict['Y'] = max(instDict['Y'],0.1) 638 instDict['SH/L'] = max(instDict['SH/L'],0.001) 639 return dataType,instDict,insVary 640 641 def GetInstParms(parmDict,Inst,varyList,Peaks): 642 instNames = Inst[3] 643 for i,name in enumerate(instNames): 644 Inst[1][i] = parmDict[name] 645 iPeak = 0 646 while True: 647 try: 648 sigName = 'sig'+str(iPeak) 649 pos = parmDict['pos'+str(iPeak)] 650 if sigName not in varyList: 651 parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W'] 652 gamName = 'gam'+str(iPeak) 653 if gamName not in varyList: 654 parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0) 655 iPeak += 1 656 except KeyError: 657 break 658 659 def InstPrint(Inst,sigDict): 593 660 print 'Instrument Parameters:' 594 661 ptfmt = "%12.6f" … … 596 663 ptstr = 'values:' 597 664 sigstr = 'esds :' 598 for i,value in enumerate(instVal): 599 ptlbls += "%s" % (insLabels[i].center(12)) 600 ptstr += ptfmt % (value) 601 if siginst[i]: 602 sigstr += ptfmt % (siginst[i]) 665 instNames = Inst[3][1:] 666 for i,name in enumerate(instNames): 667 ptlbls += "%s" % (name.center(12)) 668 ptstr += ptfmt % (Inst[1][i+1]) 669 if name in sigDict: 670 sigstr += ptfmt % (sigDict[name]) 603 671 else: 604 672 sigstr += 12*' ' … … 606 674 print ptstr 607 675 print sigstr 608 609 def peaksPrint(peaks,sigpeaks): 610 print 'Peak list=' 611 612 begin = time.time() 613 Np = len(peaks[0]) 614 DataType = inst[1][0] 615 instVal = inst[1][1:] 616 Insref = inst[2][1:] 617 insLabels = inst[3][1:] 618 Ka2 = False 619 Ioff = 3 620 if len(instVal) == 12: 621 lamratio = instVal[1]/instVal[0] 622 Ka2 = True 623 Ioff = 5 624 insref = Insref[len(Insref)7:1] #just U,V,W,X,Y,SH/L 625 for peak in peaks: 626 dip = [] 627 dip.append(tand(peak[0]/2.0)**2) 628 dip.append(tand(peak[0]/2.0)) 629 dip.append(1.0/cosd(peak[0]/2.0)) 630 dip.append(tand(peak[0]/2.0)) 631 peak.append(dip) 632 B = background[2] 633 bcof = background[3:3+B] 634 Bv = 0 635 if background[1]: 636 Bv = B 676 677 def setPeaksParms(Peaks): 678 peakNames = [] 679 peakVary = [] 680 peakVals = [] 681 names = ['pos','int','sig','gam'] 682 for i,peak in enumerate(Peaks): 683 for j in range(4): 684 peakVals.append(peak[2*j]) 685 parName = names[j]+str(i) 686 peakNames.append(parName) 687 if peak[2*j+1]: 688 peakVary.append(parName) 689 return dict(zip(peakNames,peakVals)),peakVary 690 691 def GetPeaksParms(parmDict,Peaks,varyList): 692 names = ['pos','int','sig','gam'] 693 for i,peak in enumerate(Peaks): 694 for j in range(4): 695 pos = parmDict['pos'+str(i)] 696 parName = names[j]+str(i) 697 if parName in varyList: 698 peak[2*j] = parmDict[parName] 699 elif 'sig' in parName: 700 peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W'] 701 elif 'gam' in parName: 702 peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0) 703 704 def PeaksPrint(parmDict,sigDict,varyList): 705 print 'Peak coefficients:' 706 names = ['pos','int','sig','gam'] 707 head = 15*' ' 708 for name in names: 709 head += name.center(12)+'esd'.center(12) 710 print head 711 ptfmt = {'pos':"%12.5f",'int':"%12.1f",'sig':"%12.3f",'gam':"%12.3f"} 712 for i,peak in enumerate(Peaks): 713 ptstr = ':' 714 for j in range(4): 715 name = names[j] 716 parName = name+str(i) 717 ptstr += ptfmt[name] % (parmDict[parName]) 718 if parName in varyList: 719 # ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName]) 720 ptstr += ptfmt[name] % (sigDict[parName]) 721 else: 722 # ptstr += G2IO.ValEsd(parmDict[parName],0.0) 723 ptstr += 12*' ' 724 print '%s'%(('Peak'+str(i+1)).center(8)),ptstr 725 726 def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType): 727 parmdict.update(zip(varylist,values)) 728 return np.sqrt(weights)*(ydatagetPeakProfile(parmdict,xdata,varylist,bakType)) 729 637 730 x,y,w,yc,yb,yd = data #these are numpy arrays! 638 V = [] 639 A = [] 640 swobs = 0.0 641 smin = 0.0 642 first = True 643 GoOn = True 644 Go = True 645 dlg = wx.ProgressDialog("Elapsed time","Fitting peaks to pattern",len(x), \ 646 style = wx.PD_ELAPSED_TIMEwx.PD_AUTO_HIDEwx.PD_REMAINING_TIMEwx.PD_CAN_ABORT) 647 screenSize = wx.DisplaySize() 648 Size = dlg.GetSize() 649 dlg.SetPosition(wx.Point(screenSize[0]Size[0]300,0)) 650 try: 651 i = 0 652 for xi in x : 653 Go = dlg.Update(i)[0] 654 if GoOn: 655 GoOn = Go 656 if limits[0] <= xi <= limits[1]: 657 yb[i] = 0.0 658 dp = [] 659 for j in range(B): 660 t = (xilimits[0])**j 661 yb[i] += t*bcof[j] 662 if background[1]: 663 dp.append(t) 664 yc[i] = yb[i] 665 Iv = 0 666 for j in range(6): 667 if insref[j]: 668 dp.append(0.0) 669 Iv += 1 670 for peak in peaks: 671 dip = peak[1] 672 f = pyp.pypsvfcj(peak[2],xipeak[0],peak[0],peak[4],peak[6],instVal[2],0.0) 673 yc[i] += f[0]*peak[2] 674 if f[0] > 0.0: 675 j = 0 676 if insref[0]: #U 677 dp[Bv+j] += f[3]*dip[0] 678 j += 1 679 if insref[1]: #V 680 dp[Bv+j] += f[3]*dip[1] 681 j += 1 682 if insref[2]: #W 683 dp[Bv+j] += f[3] 684 j += 1 685 if insref[3]: #X 686 dp[Bv+j] += f[4]*dip[2] 687 j += 1 688 if insref[4]: #Y 689 dp[Bv+j] += f[4]*dip[3] 690 j += 1 691 if insref[5]: #SH/L 692 dp[Bv+j] += f[5] 693 if Ka2: 694 pos2 = 2.0*asind(lamratio*sind(peak[0]/2.0)) 695 f2 = pyp.pypsvfcj(peak[2],xipos2,peak[0],peak[4],peak[6],instVal[2],0.0) 696 yc[i] += f2[0]*peak[2]*instVal[3] 697 if f[0] > 0.0: 698 j = 0 699 if insref[0]: #U 700 dp[Bv+j] += f2[3]*dip[0]*instVal[3] 701 j += 1 702 if insref[1]: #V 703 dp[Bv+j] += f2[3]*dip[1]*instVal[3] 704 j += 1 705 if insref[2]: #W 706 dp[Bv+j] += f2[3]*instVal[3] 707 j += 1 708 if insref[3]: #X 709 dp[Bv+j] += f2[4]*dip[2]*instVal[3] 710 j += 1 711 if insref[4]: #Y 712 dp[Bv+j] += f2[4]*dip[3]*instVal[3] 713 j += 1 714 if insref[5]: #SH/L 715 dp[Bv+j] += f2[5]*instVal[3] 716 for j in range(0,Np,2): 717 if peak[j+1]: dp.append(f[j/2+1]) 718 yd[i] = y[i]yc[i] 719 swobs += w[i]*y[i]**2 720 t2 = w[i]*yd[i] 721 smin += t2*yd[i] 722 if first: 723 first = False 724 M = len(dp) 725 A = np.zeros(shape=(M,M)) 726 V = np.zeros(shape=(M)) 727 A,V = pyp.buildmv(t2,w[i],M,dp,A,V) 728 i += 1 729 finally: 730 dlg.Destroy() 731 Rwp = smin/swobs 732 Rwp = math.sqrt(Rwp)*100.0 733 norm = np.diag(A) 734 for i,elm in enumerate(norm): 735 if elm <= 0.0: 736 print norm 737 return False,0,0,0,False 738 for i in xrange(len(V)): 739 norm[i] = 1.0/math.sqrt(norm[i]) 740 V[i] *= norm[i] 741 a = A[i] 742 for j in xrange(len(V)): 743 a[j] *= norm[i] 744 A = np.transpose(A) 745 for i in xrange(len(V)): 746 a = A[i] 747 for j in xrange(len(V)): 748 a[j] *= norm[i] 749 b = nl.solve(A,V) 750 A = nl.inv(A) 751 sig = np.diag(A) 752 for i in xrange(len(V)): 753 b[i] *= norm[i] 754 sig[i] *= norm[i]*norm[i] 755 sig[i] = math.sqrt(abs(sig[i])) 756 sigback = [0,0,0] 757 for j in range(Bv): 758 background[j+3] += b[j] 759 sigback.append(sig[j]) 760 backgroundPrint(background,sigback) 761 k = 0 762 delt = [] 763 if Ka2: 764 siginst = [0,0,0,0,0] 765 else: 766 siginst = [0,0,0] 767 for j in range(6): 768 if insref[j]: 769 instVal[j+Ioff] += b[Bv+k] 770 siginst.append(sig[Bv+k]) 771 delt.append(b[Bv+k]) 772 k += 1 773 else: 774 delt.append(0.0) 775 siginst.append(0.0) 776 delt.append(0.0) #dummies for azm 777 siginst.append(0.0) 778 instPrint(instVal,siginst,insLabels) 779 inst[1] = [DataType,] 780 for val in instVal: 781 inst[1].append(val) 782 B = Bv+Iv 783 for peak in peaks: 784 del peak[1] # remove dip from end 785 delsig = delt[0]*tand(peak[0]/2.0)**2+delt[1]*tand(peak[0]/2.0)+delt[2] 786 delgam = delt[3]/cosd(peak[0]/2.0)+delt[4]*tand(peak[0]/2.0) 787 for j in range(0,len(peak[:1]),2): 788 if peak[j+1]: 789 peak[j] += b[B] 790 B += 1 791 peak[4] += delsig 792 if peak[4] < 0.0: 793 print 'ERROR  negative sigma' 794 return False,0,0,0,False 795 peak[6] += delgam 796 if peak[6] < 0.0: 797 print 'ERROR  negative gamma' 798 return False,0,0,0,False 799 runtime = time.time()begin 800 data = [x,y,w,yc,yb,yd] 801 return True,smin,Rwp,runtime,GoOn 802 731 xBeg = np.searchsorted(x,Limits[0]) 732 xFin = np.searchsorted(x,Limits[1]) 733 bakType,bakDict,bakVary = SetBackgroundParms(Background) 734 dataType,insDict,insVary = SetInstParms(Inst) 735 peakDict,peakVary = setPeaksParms(Peaks) 736 parmDict = {} 737 parmDict.update(bakDict) 738 parmDict.update(insDict) 739 parmDict.update(peakDict) 740 varyList = bakVary+insVary+peakVary 741 while True: 742 begin = time.time() 743 values = np.array(ValuesOut(parmDict, varyList)) 744 if FitPgm == 'LSQ': 745 result = so.leastsq(errPeakProfile,values, 746 args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType),full_output=True) 747 runtime = time.time()begin 748 print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFinxBeg,' Number of parameters: ',len(varyList) 749 print "%s%8.3f%s " % ('fitpeak time =',runtime,'s') 750 ValuesIn(parmDict, varyList, result[0]) 751 chisq = np.sum(errPeakProfile(result[0],x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType)**2) 752 Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100. #to % 753 GOF = chisq/(xFinxBeglen(varyList)) 754 print "%s%7.2f%s%12.6g%s%6.2f" % ('Rwp = ',Rwp,'%, chi**2 = ',chisq,' reduced chi**2 = ',GOF) 755 try: 756 sig = np.sqrt(np.diag(result[1])*chisq) 757 break #refinement succeeded  finish up! 758 except ValueError: #result[1] is None on singular matrix 759 print 'Refinement failed  singular matrix' 760 Ipvt = result[2]['ipvt'] 761 for i,ipvt in enumerate(Ipvt): 762 if not np.sum(result[2]['fjac'],axis=1)[i]: 763 print 'Removing parameter: ',varyList[ipvt1] 764 del(varyList[ipvt1]) 765 break 766 elif FitPgm == 'BFGS': 767 print 'Other program here' 768 return 769 770 sigDict = dict(zip(varyList,sig)) 771 yb[xBeg:xFin] = getBackground(parmDict,bakType,x[xBeg:xFin]) 772 yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType) 773 yd[xBeg:xFin] = y[xBeg:xFin]yc[xBeg:xFin] 774 GetBackgroundParms(parmDict,Background) 775 BackgroundPrint(Background,sigDict) 776 GetInstParms(parmDict,Inst,varyList,Peaks) 777 InstPrint(Inst,sigDict) 778 GetPeaksParms(parmDict,Peaks,varyList) 779 PeaksPrint(parmDict,sigDict,varyList) 780 803 781 def CalcPDF(data,inst,xydata): 804 782 auxPlot = [] … … 881 859 return auxPlot 882 860 861 #testing data 862 import plot 863 NeedTestData = True 864 def TestData(): 865 # global NeedTestData 866 NeedTestData = False 867 global bakType 868 bakType = 'chebyschev' 869 global xdata 870 xdata = np.linspace(4.0,40.0,36000) 871 global parmDict0 872 parmDict0 = { 873 'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0, 874 'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0, 875 'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0, 876 'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0, 877 'U':1.163,'V':0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002, 878 'Back0':5.384,'Back1':0.015,'Back2':.004, 879 } 880 global parmDict1 881 parmDict1 = { 882 'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0, 883 'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0, 884 'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0, 885 'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0, 886 'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0, 887 'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0, 888 'U':22.75,'V':17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002, 889 'Back0':36.897,'Back1':0.508,'Back2':.006, 890 'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5, 891 } 892 global varyList 893 varyList = [] 894 895 def test0(): 896 if NeedTestData: TestData() 897 msg = 'test ' 898 gplot = plotter.add('FCJVoigt, 11BM').gca() 899 gplot.plot(xdata,getBackground(parmDict0,bakType,xdata)) 900 gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType)) 901 fplot = plotter.add('FCJVoigt, Ka1+2').gca() 902 fplot.plot(xdata,getBackground(parmDict1,bakType,xdata)) 903 fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType)) 904 time0 = time.time() 905 for i in range(100): 906 y = getPeakProfile(parmDict1,xdata,varyList,bakType) 907 print '100+6*Ka12 peaks=3600 peaks',time.time()time0 908 909 910 if __name__ == '__main__': 911 global plotter 912 plotter = plot.PlotNotebook() 913 test0() 914 print "OK" 915 plotter.StartEventLoop()
Note: See TracChangeset
for help on using the changeset viewer.