Changeset 1252
 Timestamp:
 Mar 18, 2014 1:45:39 PM (9 years ago)
 Location:
 trunk
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIplot.py
r1250 r1252 890 890 D = xye[5]Ymax*G2frame.delOffset 891 891 elif 'SASD' in plottype: 892 D = xye[4] Ymax*G2frame.delOffset892 D = xye[4] 893 893 Plot.set_yscale("log",nonposy='mask') 894 894 Plot.set_ylim(bottom=np.min(np.trim_zeros(Y))/2.,top=np.max(Y)*2.) … … 907 907 else: 908 908 Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False) 909 Plot.plot(X,D,colors[(N+1)%6],picker=False) 909 910 Plot.plot(X,Z,colors[(N+1)%6],picker=False) 910 911 elif G2frame.Weight and 'PWDR' in plottype: 
trunk/GSASIIpwdGUI.py
r1250 r1252 86 86 'logBins':True,'Method':'MaxEnt','Distribution':[], 87 87 'Shape':['Spheroid',1.0],'MaxEnt':{'Niter':100,'Precision':0.01,'Sky':3}, 88 'IPG':{'Niter':100,'Approach':0.8 },'Reg':{},},88 'IPG':{'Niter':100,'Approach':0.8,'Power':1},'Reg':{},}, 89 89 'Unified':{'Levels':[],}, 90 90 'Particle':{'Levels':[],}, … … 2475 2475 if isinstance(data['Size']['MaxEnt']['Sky'],float): 2476 2476 data['Size']['MaxEnt']['Sky'] = 3 2477 if 'Power' not in data['Size']['IPG']: 2478 data['Size']['IPG']['Power'] = 1 2477 2479 #end patches 2478 2480 … … 2483 2485 def OnFitModel(event): 2484 2486 print 'fit model for '+data['Current'] 2487 if not any(Sample['Contrast']): 2488 G2frame.ErrorDialog('No contrast; your sample is a vacuum!', 2489 'You need to define a scattering substance!\n'+ \ 2490 ' Do Substances and then Sample parameters') 2491 return 2485 2492 if data['Current'] == 'Size dist.': 2486 2493 G2sasd.SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data) … … 2501 2508 Obj.SetValue(fmt%(value)) 2502 2509 data[itemKey][ind] = value 2510 if itemKey == 'Back': 2511 Profile[4][:] = value 2512 G2plt.PlotPatterns(G2frame,plotType='SASD',newPlot=True) 2503 2513 2504 2514 def OnCheckBox(event): … … 2587 2597 sizeSizer.Add((5,5),0) 2588 2598 fitSizer = wx.BoxSizer(wx.HORIZONTAL) 2589 methods = ['MaxEnt', 'IPG',]2599 methods = ['MaxEnt',] #'IPG',] 2590 2600 fitSizer.Add(wx.StaticText(G2frame.dataDisplay,label='Fitting method: '),0,WACV) 2591 2601 method = wx.ComboBox(G2frame.dataDisplay,value=data['Size']['Method'],choices=methods, … … 2609 2619 floor.Bind(wx.EVT_COMBOBOX,OnIntVal) 2610 2620 fitSizer.Add(floor,0,WACV) 2621 elif 'IPG' in data['Size']['Method']: 2622 fitSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Q power weight (1 for sigma): '),0,WACV) 2623 choices = ['1','0','1','2','3','4'] 2624 power = wx.ComboBox(G2frame.dataDisplay,value=str(data['Size']['IPG']['Power']),choices=choices, 2625 style=wx.CB_READONLYwx.CB_DROPDOWN) 2626 Indx[power.GetId()] = [data['Size']['IPG'],'Power',2] 2627 power.Bind(wx.EVT_COMBOBOX,OnIntVal) 2628 fitSizer.Add(power,0,WACV) 2611 2629 sizeSizer.Add(fitSizer,0) 2612 2630 … … 2668 2686 G2gd.HorizontalLine(mainSizer,G2frame.dataDisplay) 2669 2687 if 'Size' in data['Current']: 2688 if 'MaxEnt' in data['Size']['Method']: 2689 Status.SetStatusText('Size distribution by Maximum entropy') 2690 elif 'IPG' in data['Size']['Method']: 2691 Status.SetStatusText('Size distribution by InteriorPoint Gradient') 2670 2692 mainSizer.Add(SizeSizer()) 2671 2693 elif 'Particle' in data['Current']: … … 2681 2703 backVal.Bind(wx.EVT_KILL_FOCUS,OnValueChange) 2682 2704 backSizer.Add(backVal,0,WACV) 2683 backVar = wx.CheckBox(G2frame.dataDisplay,label=' Apply?')2705 backVar = wx.CheckBox(G2frame.dataDisplay,label='Refine?') 2684 2706 Indx[backVar.GetId()] = [data['Back'],1] 2685 2707 backVar.SetValue(data['Back'][1]) 
trunk/GSASIIsasd.py
r1250 r1252 63 63 def SphereFF(Q,R,args=()): 64 64 ''' Compute hard sphere form factor  can use numpy arrays 65 param float:Q Q value array (usually in A1) 66 param float:R sphere radius (Usually in A  must match Q1 units) 65 param float Q: Q value array (usually in A1) 66 param float R: sphere radius (Usually in A  must match Q1 units) 67 param array args: ignored 67 68 returns float: form factors as array as needed 68 69 ''' … … 73 74 ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid) 74 75  can use numpy arrays for R & AR; will return corresponding numpy array 75 param float :QQ value array (usually in A1)76 param float Q : Q value array (usually in A1) 76 77 param float R: radius along 2 axes of spheroid 77 param float AR: aspect ratio so 3rd axis = R*AR78 param array args: [float AR]: aspect ratio so 3rd axis = R*AR 78 79 returns float: form factors as array as needed 79 80 ''' … … 89 90 def CylinderFF(Q,R,args): 90 91 ''' Compute form factor for cylinders  can use numpy arrays 91 param float : QQ value array (A1)92 param float : Rcylinder radius (A)93 param float: Lcylinder length (A)92 param float Q: Q value array (A1) 93 param float R: cylinder radius (A) 94 param array args: [float L]: cylinder length (A) 94 95 returns float: form factor 95 96 ''' … … 107 108 def CylinderDFF(Q,L,args): 108 109 ''' Compute form factor for cylinders  can use numpy arrays 109 param float : QQ value array (A1)110 param float : Lcylinder half length (A)111 param float: R cylinder diameter(A)110 param float Q: Q value array (A1) 111 param float L: cylinder half length (A) 112 param array args: [float R]: cylinder radius (A) 112 113 returns float: form factor 113 114 ''' … … 117 118 def CylinderARFF(Q,R,args): 118 119 ''' Compute form factor for cylinders  can use numpy arrays 119 param float : QQ value array (A1)120 param float : Rcylinder radius (A)121 param float: ARcylinder aspect ratio = L/D = L/2R120 param float Q: Q value array (A1) 121 param float R: cylinder radius (A) 122 param array args: [float AR]: cylinder aspect ratio = L/D = L/2R 122 123 returns float: form factor 123 124 ''' … … 126 127 127 128 def UniSphereFF(Q,R,args=0): 129 ''' Compute form factor for unified sphere  can use numpy arrays 130 param float Q: Q value array (A1) 131 param float R: cylinder radius (A) 132 param array args: ignored 133 returns float: form factor 134 ''' 128 135 Rg = np.sqrt(3./5.)*R 129 B = 1.62/(Rg**4) #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?136 B = np.pi*1.62/(Rg**4) #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense? 130 137 QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg/np.sqrt(6)))**3 131 138 return np.sqrt(np.exp((Q[:,np.newaxis]**2*Rg**2)/3.)+(B/QstV**4)) 132 139 133 140 def UniRodFF(Q,R,args): 141 ''' Compute form factor for unified rod  can use numpy arrays 142 param float Q: Q value array (A1) 143 param float R: cylinder radius (A) 144 param array args: [float R]: cylinder radius (A) 145 returns float: form factor 146 ''' 134 147 L = args[0] 135 148 Rg2 = np.sqrt(R**2/2+L**2/12) … … 145 158 146 159 def UniRodARFF(Q,R,args): 160 ''' Compute form factor for unified rod of fixed aspect ratio  can use numpy arrays 161 param float Q: Q value array (A1) 162 param float R: cylinder radius (A) 163 param array args: [float AR]: cylinder aspect ratio = L/D = L/2R 164 returns float: form factor 165 ''' 147 166 AR = args[0] 148 167 return UniRodFF(Q,R,[AR*R,]) 149 168 150 169 def UniDiskFF(Q,R,args): 170 ''' Compute form factor for unified disk  can use numpy arrays 171 param float Q: Q value array (A1) 172 param float R: cylinder radius (A) 173 param array args: [float T]: disk thickness (A) 174 returns float: form factor 175 ''' 151 176 T = args[0] 152 177 Rg2 = np.sqrt(R**2/2.+T**2/12.) … … 163 188 164 189 def UniTubeFF(Q,R,args): 190 ''' Compute form factor for unified disk  can use numpy arrays 191 param float Q: Q value array (A1) 192 param float R: cylinder radius (A) 193 param array args: [float L,T]: tube length & wall thickness(A) 194 returns float: form factor 195 ''' 165 196 L,T = args[:2] 166 197 Ri = RT … … 186 217 ''' Compute volume of sphere 187 218  numpy array friendly 188 param float:R sphere radius 219 param float R: sphere radius 220 param array args: ignored 189 221 returns float: volume 190 222 ''' … … 195 227  numpy array friendly 196 228 param float R: radius along 2 axes of spheroid 197 param float AR: aspect ratio so radius of 3rd axis = R*AR229 param array args: [float AR]: aspect ratio so radius of 3rd axis = R*AR 198 230 returns float: volume 199 231 ''' … … 204 236 ''' Compute cylinder volume for radius & length 205 237  numpy array friendly 206 param float : Rdiameter (A)207 param float: Llength (A)238 param float R: diameter (A) 239 param array args: [float L]: length (A) 208 240 returns float:volume (A^3) 209 241 ''' … … 215 247  numpy array friendly 216 248 param float: L half length (A) 217 param float: Ddiameter (A)249 param array args: [float D]: diameter (A) 218 250 returns float:volume (A^3) 219 251 ''' … … 225 257  numpy array friendly 226 258 param float: R radius (A) 227 param float: AR=L/D=L/2R aspect ratio259 param array args: [float AR]: =L/D=L/2R aspect ratio 228 260 returns float:volume 229 261 ''' … … 234 266 ''' Compute volume of sphere 235 267  numpy array friendly 236 param float:R sphere radius 268 param float R: sphere radius 269 param array args: ignored 237 270 returns float: volume 238 271 ''' … … 242 275 ''' Compute cylinder volume for radius & length 243 276  numpy array friendly 244 param float : Rdiameter (A)245 param float: Llength (A)277 param float R: diameter (A) 278 param array args: [float L]: length (A) 246 279 returns float:volume (A^3) 247 280 ''' … … 250 283 251 284 def UniRodARVol(R,args): 285 ''' Compute rod volume for radius & aspect ratio 286  numpy array friendly 287 param float R: diameter (A) 288 param array args: [float AR]: =L/D=L/2R aspect ratio 289 returns float:volume (A^3) 290 ''' 252 291 AR = args[0] 253 292 return CylinderARVol(R,[AR,]) 254 293 255 294 def UniDiskVol(R,args): 295 ''' Compute disk volume for radius & thickness 296  numpy array friendly 297 param float R: diameter (A) 298 param array args: [float T]: thickness 299 returns float:volume (A^3) 300 ''' 256 301 T = args[0] 257 302 return CylinderVol(R,[T,]) … … 260 305 ''' Compute tube volume for radius, length & wall thickness 261 306  numpy array friendly 262 param float: R diameter (A) 263 param float: L length (A) 264 param float: T tube wall thickness (A) 307 param float R: diameter (A) 308 param array args: [float L,T]: tube length & wall thickness(A) 265 309 returns float: volume (A^3) of tube wall 266 310 ''' … … 317 361 pass 318 362 319 def MaxEnt_SB(datum, sigma, base, IterMax, G, image_to_data=None, data_to_image=None, report=False):363 def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False): 320 364 ''' 321 365 do the complete Maximum Entropy algorithm of Skilling and Bryan … … 323 367 :param float datum[]: 324 368 :param float sigma[]: 369 :param float[][] G: transformation matrix 325 370 :param float base[]: 326 371 :param int IterMax: 327 :param float[][] G: transformation matrix328 372 :param obj image_to_data: opus function (defaults to opus) 329 373 :param obj data_to_image: tropus function (defaults to tropus) … … 617 661 return chisq,f,image_to_data(f, G) # no solution after IterMax iterations 618 662 663 664 ############################################################################### 665 #### IPG/TNNLS Routines 666 ############################################################################### 667 668 def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=1,report=False): 669 ''' An implementation of the InteriorPoint Gradient method of 670 Michael Merritt & Yin Zhang, Technical Report TR0408, Dept. of Comp. and 671 Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at 672 http://www.caam.rice.edu/caam/trs/2004/TR0408.pdf 673 Problem addressed: Total NonNegative Least Squares (TNNLS) 674 :param float datum[]: 675 :param float sigma[]: 676 :param float[][] G: transformation matrix 677 :param int IterMax: 678 :param float Qvec: data positions for Power = 04 679 :param float approach: 0.8 default fitting parameter 680 :param int Power: 04 for Q^Power weighting, 1 to use input sigma 681 682 ''' 683 if Power < 0: 684 GmatE = G/sigma[:np.newaxis] 685 IntE = datum/sigma 686 pwr = 0 687 QvecP = np.ones_like(datum) 688 else: 689 GmatE = G[:] 690 IntE = datum[:] 691 pwr = Power 692 QvecP = Qvec**pwr 693 Amat = GmatE*QvecP[:np.newaxis] 694 AAmat = np.inner(Amat,Amat) 695 Bvec = datum*QvecP 696 Xw = np.ones_like(Bins)*1.e6 697 calc = np.dot(G.T,Xw) 698 nIter = 0 699 err = 10. 700 while (nIter<IterMax) and (err > 1.): 701 #Step 1 in M&Z paper: 702 Qk = np.dot(AAmat,Xw)np.dot(Amat,Bvec) 703 Dk = Xw/np.dot(AAmat,Xw) 704 Pk = Dk*Qk 705 #Step 2 in M&Z paper: 706 alpSt = np.dot(Pk,Qk)/np.dot(Pk,np.dot(AAmat,Pk)) 707 alpWv = Xw/Pk 708 AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),Pk[i]) for i in range(Pk.shape[0])] 709 #Step 3 in M&Z paper: 710 shift = AkSt*Pk 711 print np.sum(shift**2) 712 Xw += shift 713 #done IPG; now results 714 nIter += 1 715 calc = np.dot(G.T,Xw) 716 chisq = np.sum(((datumcalc)/sigma)**2) 717 err = chisq/len(datum) 718 if report: 719 print ' Iteration: %d, chisq: %.3g'%(nIter,chisq) 720 return chisq,Xw,calc 721 722 ############################################################################### 723 #### SASD Utilities 724 ############################################################################### 725 726 def SetScale(Data,refData): 727 Profile,Limits,Sample = Data 728 refProfile,refLimits,refSample = refData 729 x,y = Profile[:2] 730 rx,ry = refProfile[:2] 731 Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]]) 732 Fin = np.min([rx[1],x[1],Limits[1][1],refLimits[1][1]]) 733 iBeg = np.searchsorted(x,Beg) 734 iFin = np.searchsorted(x,Fin) 735 sum = np.sum(y[iBeg:iFin]) 736 refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0)) 737 Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum 738 739 ############################################################################### 740 #### Size distribution 741 ############################################################################### 742 743 def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data): 744 shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol], 745 'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol], 746 'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol], 747 'Unified disk':[UniDiskFF,UniDiskVol]} 748 Shape = data['Size']['Shape'][0] 749 Parms = data['Size']['Shape'][1:] 750 if data['Size']['logBins']: 751 Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']), 752 data['Size']['Nbins']+1,True)/2. #make radii 753 else: 754 Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'], 755 data['Size']['Nbins']+1,True)/2. #make radii 756 Dbins = np.diff(Bins) 757 Bins = Bins[:1]+Dbins/2. 758 Contrast = Sample['Contrast'][1] 759 Scale = Sample['Scale'][0] 760 Sky = 10**data['Size']['MaxEnt']['Sky'] 761 BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast #How about *Scale/Contrast? 762 Back = data['Back'] 763 Q,Io,wt,Ic,Ib = Profile 764 Qmin = Limits[1][0] 765 Qmax = Limits[1][1] 766 wtFactor = ProfDict['wtFactor'] 767 Ibeg = np.searchsorted(Q,Qmin) 768 Ifin = np.searchsorted(Q,Qmax) 769 BinMag = np.zeros_like(Bins) 770 Ic[:] = 0. 771 Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms) 772 if 'MaxEnt' == data['Size']['Method']: 773 chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]Back[0], 774 Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack, 775 data['Size']['MaxEnt']['Niter'],report=True) 776 elif 'IPG' == data['Size']['Method']: 777 chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]), 778 Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8, 779 Power=data['Size']['IPG']['Power'],report=True) 780 Ib[:] = Back[0] 781 Ic[Ibeg:Ifin] += Back[0] 782 print ' Final chi^2: %.3f'%(chisq) 783 Vols = shapes[Shape][1](Bins,Parms) 784 data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)] 785 619 786 620 787 ################################################################################ … … 685 852 if __name__ == '__main__': 686 853 tests() 687 688 ############################################################################### 689 #### SASD Utilities 690 ############################################################################### 691 692 def SetScale(Data,refData): 693 Profile,Limits,Sample = Data 694 refProfile,refLimits,refSample = refData 695 x,y = Profile[:2] 696 rx,ry = refProfile[:2] 697 Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]]) 698 Fin = np.min([rx[1],x[1],Limits[1][1],refLimits[1][1]]) 699 iBeg = np.searchsorted(x,Beg) 700 iFin = np.searchsorted(x,Fin) 701 sum = np.sum(y[iBeg:iFin]) 702 refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0)) 703 Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum 704 705 ############################################################################### 706 #### Size distribution 707 ############################################################################### 708 709 def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data): 710 shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol], 711 'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol], 712 'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol], 713 'Unified disk':[UniDiskFF,UniDiskVol]} 714 Shape = data['Size']['Shape'][0] 715 Parms = data['Size']['Shape'][1:] 716 if data['Size']['logBins']: 717 Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']), 718 data['Size']['Nbins']+1,True)/2. #make radii 719 else: 720 Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'], 721 data['Size']['Nbins']+1,True)/2. #make radii 722 Dbins = np.diff(Bins) 723 Bins = Bins[:1]+Dbins/2. 724 Contrast = Sample['Contrast'][1] 725 Scale = Sample['Scale'][0] 726 Sky = 10**data['Size']['MaxEnt']['Sky'] 727 BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast #How about *Scale/Contrast? 728 Back = data['Back'] 729 Q,Io,wt,Ic,Ib = Profile[:5] 730 Qmin = Limits[1][0] 731 Qmax = Limits[1][1] 732 wtFactor = ProfDict['wtFactor'] 733 Ibeg = np.searchsorted(Q,Qmin) 734 Ifin = np.searchsorted(Q,Qmax) 735 Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms) 736 chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]Back[0], 737 Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),BinsBack, 738 data['Size']['MaxEnt']['Niter'],Gmat,report=True) 739 if Back[1]: 740 Ib = Back[0] 741 Ic[Ibeg:Ifin] += Back[0] 742 print ' Final chi^2: %.3f'%(chisq) 743 Vols = shapes[Shape][1](Bins,Parms) 744 data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)] 745 746 854
Note: See TracChangeset
for help on using the changeset viewer.