Changeset 2757
- Timestamp:
- Mar 21, 2017 11:01:51 AM (6 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIpwd.py
r2755 r2757 1913 1913 ################################################################################ 1914 1914 1915 def REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data): 1915 def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data): 1916 print 'fit REFD data' 1917 1918 def GetModelParms(): 1919 parmDict = {} 1920 varyList = [] 1921 values = [] 1922 for parm in ['Scale','FltBack']: 1923 parmDict[parm] = data[parm][0] 1924 if data[parm][1]: 1925 varyList.append(parm) 1926 values.append(data[parm][0]) 1927 parmDict['nLayers'] = len(data['Layers']) 1928 for ilay,layer in enumerate(data['Layers']): 1929 name = layer['Name'] 1930 cid = str(ilay)+';' 1931 for parm in ['Thick','Rough','DenMul']: 1932 parmDict[cid+parm] = layer.get(parm,[0.,False])[0] 1933 if layer.get(parm,[0.,False])[1]: 1934 varyList.append(cid+parm) 1935 values.append(layer[parm][0]) 1936 parmDict[cid+'rho'] = Substances[name]['Scatt density'] 1937 parmDict[cid+'irho'] = Substances[name].get('XImag density',0.) 1938 return parmDict,varyList,values 1939 1940 def SetModelParms(): 1941 line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale']) 1942 if 'Scale' in varyList: 1943 data['Scale'] = parmDict['Scale'] 1944 line += ' esd: %.4g'%(sigDict['Scale']) 1945 print line 1946 line = ' Flat background: %15.4g'%(parmDict['FltBack']) 1947 if 'FltBack' in varyList: 1948 data['FltBack'] = parmDict['FltBack'] 1949 line += ' esd: %15.3g'%(sigDict['FltBack']) 1950 print line 1951 for ilay,layer in enumerate(data['Layers']): 1952 name = layer['Name'] 1953 print ' Parameters for layer: %d %s'%(ilay,name) 1954 cid = str(ilay)+';' 1955 line = ' ' 1956 line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul']) 1957 line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)**parmDict[cid+'DenMul']) 1958 for parm in ['Thick','Rough','DenMul']: 1959 if parm in layer: 1960 layer[parm][0] = parmDict[cid+parm] 1961 line += ' %s: %.3f'%(parm,layer[parm][0]) 1962 if cid+parm in varyList: 1963 line += ' esd: %.3g'%(sigDict[cid+parm]) 1964 print line 1965 print line2 1966 1967 def calcREFD(values,Q,Io,wt,parmDict,varyList): 1968 parmDict.update(zip(varyList,values)) 1969 M = np.sqrt(wt)*(getREFD(Q,parmDict)-Io) 1970 return M 1971 1972 def getREFD(Q,parmDict): 1973 Ic = np.ones_like(Q)*parmDict['FltBack'] 1974 Scale = parmDict['Scale'] 1975 Nlayers = parmDict['nLayers'] 1976 depth = np.zeros(Nlayers) 1977 rho = np.zeros(Nlayers) 1978 irho = np.zeros(Nlayers) 1979 sigma = np.zeros(Nlayers) 1980 for ilay in range(Nlayers): 1981 cid = str(ilay)+';' 1982 depth[ilay] = parmDict[cid+'Thick'] 1983 sigma[ilay] = parmDict[cid+'Rough'] 1984 rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul'] 1985 irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul'] 1986 A,B = abeles(0.5*Q,depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles 1987 Ic += (A**2+B**2)*Scale 1988 return Ic 1989 1990 Q,Io,wt,Ic,Ib,Ifb = Profile[:6] 1991 Qmin = Limits[1][0] 1992 Qmax = Limits[1][1] 1993 wtFactor = ProfDict['wtFactor'] 1994 Ibeg = np.searchsorted(Q,Qmin) 1995 Ifin = np.searchsorted(Q,Qmax)+1 #include last point 1996 Ic[:] = 0 1997 parmDict,varyList,values = GetModelParms() 1998 if varyList: 1999 result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8, #ftol=Ftol, 2000 args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)) 2001 parmDict.update(zip(varyList,result[0])) 2002 chisq = np.sum(result[2]['fvec']**2) 2003 ncalc = result[2]['nfev'] 2004 covM = result[1] 2005 else: #nothing varied 2006 M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList) 2007 chisq = np.sum(M**2) 2008 ncalc = 0 2009 covM = [] 2010 sig = [] 2011 sigDict = {} 2012 result = [] 2013 Rvals = {} 2014 Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100. #to % 2015 Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList)) #reduced chi^2 2016 Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],parmDict) 2017 Ib[Ibeg:Ifin] = parmDict['FltBack'] 2018 Msg = 'Failed to converge' 2019 try: 2020 Nans = np.isnan(result[0]) 2021 if np.any(Nans): 2022 name = varyList[Nans.nonzero(True)[0]] 2023 Msg += ' Nan result for '+name+'!' 2024 raise ValueError 2025 Negs = np.less_equal(result[0],0.) 2026 if np.any(Negs): 2027 indx = Negs.nonzero() 2028 name = varyList[indx[0][0]] 2029 if name != 'FltBack': 2030 Msg += ' negative coefficient for '+name+'!' 2031 raise ValueError 2032 if len(covM): 2033 sig = np.sqrt(np.diag(covM)*Rvals['GOF']) 2034 sigDict = dict(zip(varyList,sig)) 2035 print ' Results of reflectometry data modelling fit:' 2036 print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList) 2037 print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']) 2038 SetModelParms() 2039 covMatrix = covM*Rvals['GOF'] 2040 return True,result,varyList,sig,Rvals,covMatrix,parmDict,'' 2041 except (ValueError,TypeError): #when bad LS refinement; covM missing or with nans 2042 print Msg 2043 return False,0,0,0,0,0,0,Msg 2044 2045 def REFDModelFxn(Profile,Inst,Limits,Substances,data): 1916 2046 1917 2047 Q,Io,wt,Ic,Ib,Ifb = Profile[:6] -
trunk/GSASIIpwdGUI.py
r2755 r2757 4706 4706 def OnFitModel(event): 4707 4707 4708 print 'fit model' 4709 # SaveState() 4710 G2pwd.REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data) 4711 G2plt.PlotPatterns(G2frame,plotType='REFD') 4712 event.Skip() 4708 SaveState() 4709 G2pwd.REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data) 4710 wx.CallAfter(G2plt.PlotPatterns,G2frame,plotType='REFD') 4711 wx.CallAfter(UpdateREFDModelsGrid,G2frame,data) 4713 4712 4714 4713 def OnFitModelAll(event): … … 4720 4719 data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame, 4721 4720 G2frame.PatternId,'Models')) 4722 G2frame.dataFrame.R efdUndo.Enable(False)4721 G2frame.dataFrame.REFDUndo.Enable(False) 4723 4722 UpdateREFDModelsGrid(G2frame,data) 4724 G2pwd.REFDModelFxn(Profile, ProfDict,Inst,Limits,Substances,data)4723 G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data) 4725 4724 4726 4725 def DoUnDo(): 4727 4726 print 'Undo last refinement' 4728 file = open(G2frame.undo sasd,'rb')4727 file = open(G2frame.undorefd,'rb') 4729 4728 PatternId = G2frame.PatternId 4730 4729 G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Models'),cPickle.load(file)) 4731 print ' Model srecovered'4730 print ' Model recovered' 4732 4731 file.close() 4733 4732 … … 4736 4735 file = open(G2frame.undorefd,'wb') 4737 4736 PatternId = G2frame.PatternId 4738 for item in ['Models']: 4739 cPickle.dump(G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,item)),file,1) 4737 cPickle.dump(G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Models')),file,1) 4740 4738 file.close() 4741 G2frame.dataFrame.R efdUndo.Enable(True)4739 G2frame.dataFrame.REFDUndo.Enable(True) 4742 4740 4743 4741 def ControlSizer(): … … 4785 4783 if invalid: 4786 4784 return 4787 G2pwd.REFDModelFxn(Profile, ProfDict,Inst,Limits,Substances,data)4785 G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data) 4788 4786 G2plt.PlotPatterns(G2frame,plotType='REFD') 4789 4787 … … 4815 4813 data['Layers'][item]['Rough'] = [0.,False] 4816 4814 data['Layers'][item]['Thick'] = [1.,False] 4817 G2pwd.REFDModelFxn(Profile, ProfDict,Inst,Limits,Substances,data)4815 G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data) 4818 4816 G2plt.PlotPatterns(G2frame,plotType='REFD') 4819 4817 wx.CallAfter(UpdateREFDModelsGrid,G2frame,data) … … 4828 4826 ind = Indx[Obj.GetId()] 4829 4827 data['Layers'].insert(ind+1,{'Name':'vacuum','DenMul':[1.0,False],}) 4830 G2pwd.REFDModelFxn(Profile, ProfDict,Inst,Limits,Substances,data)4828 G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data) 4831 4829 G2plt.PlotPatterns(G2frame,plotType='REFD') 4832 4830 wx.CallAfter(UpdateREFDModelsGrid,G2frame,data) … … 4836 4834 ind = Indx[Obj.GetId()] 4837 4835 del data['Layers'][ind] 4838 G2pwd.REFDModelFxn(Profile, ProfDict,Inst,Limits,Substances,data)4836 G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data) 4839 4837 G2plt.PlotPatterns(G2frame,plotType='REFD') 4840 4838 wx.CallAfter(UpdateREFDModelsGrid,G2frame,data) … … 4843 4841 if invalid: 4844 4842 return 4845 G2pwd.REFDModelFxn(Profile, ProfDict,Inst,Limits,Substances,data)4843 G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data) 4846 4844 G2plt.PlotPatterns(G2frame,plotType='REFD') 4847 4845 wx.CallAfter(UpdateREFDModelsGrid,G2frame,data) … … 5508 5506 data['delt-G(R)'][2] += ('-\n'+subData[2]) 5509 5507 G2plt.PlotISFG(G2frame,data,newPlot=True,plotType='delt-G(R)') 5508 wx.CallAfter(UpdatePDFGrid,G2frame,data) 5509 5510 def OnMult(invalid,value,tc): 5511 if invalid: return 5512 id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,data['diffGRname']) 5513 pId = G2gd.GetPatternTreeItemId(G2frame,id,'PDF Controls') 5514 subData = G2frame.PatternTree.GetItemPyData(pId)['G(R)'] 5515 data['delt-G(R)'][1] = np.array([subData[1][0],data['G(R)'][1][1]-data['diffMult']*subData[1][1]]) 5516 G2plt.PlotISFG(G2frame,data,newPlot=True,plotType='delt-G(R)') 5510 5517 5511 5518 diffSizer = wx.BoxSizer(wx.HORIZONTAL) … … 5515 5522 style=wx.CB_READONLY|wx.CB_DROPDOWN) 5516 5523 grName.Bind(wx.EVT_COMBOBOX,OnSelectGR) 5517 diffSizer.Add(grName,0,) 5524 diffSizer.Add(grName,0,WACV) 5525 if data['diffGRname']: 5526 diffSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Mult: '),0,WACV) 5527 mult = G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'diffMult',nDig=(10,3), 5528 typeHint=float,OnLeave=OnMult) 5529 diffSizer.Add(mult,0,WACV) 5530 OnMult(False,None,None) 5518 5531 return diffSizer 5519 5532 … … 5760 5773 if 'diffGRname' not in data: 5761 5774 data['diffGRname'] = '' 5775 if 'diffMult' not in data: 5776 data['diffMult'] = 1.0 5762 5777 if G2frame.dataDisplay: 5763 5778 G2frame.dataFrame.Clear()
Note: See TracChangeset
for help on using the changeset viewer.