Changeset 1790

Timestamp:

Apr 15, 2015 4:20:11 PM (10 years ago)

Author:

vondreele

Message:

remove unused FORPI from G2strMath
New fitTexture routine to use seq refinement results - seem ok, but needs checking
bring DefaultControls? up to date in G2obj

Location:

trunk

Files:

• GSASIImath.py (modified) (2 diffs)
• GSASIIobj.py (modified) (1 diff)
• GSASIIphsGUI.py (modified) (3 diffs)
• GSASIIstrMath.py (modified) (2 diffs)

trunk/GSASIImath.py ¶

@@ -33 +33 @@
 import GSASIIpwd as G2pwd
 import numpy.fft as fft
+import scipy.optimize as so
 import pypowder as pyd
 
@@ -1592 +1593 @@
 ################################################################################
 
-def FitTexture(General,SamAngs,refData):
-    print 'fit texture for '+General["Name"]
+def FitTexture(General,Gangls,refData):
+    
+    def printSpHarm(textureData,SHtextureSig):
+        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
+        names = ['omega','chi','phi']
+        namstr = '  names :'
+        ptstr =  '  values:'
+        sigstr = '  esds  :'
+        for name in names:
+            namstr += '%12s'%('Sample '+name)
+            ptstr += '%12.3f'%(textureData['Sample '+name][1])
+            if 'Sample '+name in SHtextureSig:
+                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
+            else:
+                sigstr += 12*' '
+        print namstr
+        print ptstr
+        print sigstr
+        print '\n Texture coefficients:'
+        SHcoeff = textureData['SH Coeff'][1]
+        SHkeys = SHcoeff.keys()
+        nCoeff = len(SHcoeff)
+        nBlock = nCoeff/10+1
+        iBeg = 0
+        iFin = min(iBeg+10,nCoeff)
+        for block in range(nBlock):
+            namstr = '  names :'
+            ptstr =  '  values:'
+            sigstr = '  esds  :'
+            for name in SHkeys[iBeg:iFin]:
+                if 'C' in name:
+                    namstr += '%12s'%(name)
+                    ptstr += '%12.3f'%(SHcoeff[name])
+                    if name in SHtextureSig:
+                        sigstr += '%12.3f'%(SHtextureSig[name])
+                    else:
+                        sigstr += 12*' '
+            print namstr
+            print ptstr
+            print sigstr
+            iBeg += 10
+            iFin = min(iBeg+10,nCoeff)
+            
+    def Dict2Values(parmdict, varylist):
+        '''Use before call to leastsq to setup list of values for the parameters 
+        in parmdict, as selected by key in varylist'''
+        return [parmdict[key] for key in varylist] 
+        
+    def Values2Dict(parmdict, varylist, values):
+        ''' Use after call to leastsq to update the parameter dictionary with 
+        values corresponding to keys in varylist'''
+        parmdict.update(zip(varylist,values))
+        
+    def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList):
+        parmDict.update(zip(varyList,values))
+        Mat = np.empty(0)
+        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
+        for hist in Gangls.keys():
+            mat = []
+            for ref in refData[hist]:
+#                wt = np.sqrt(ref[4])
+                wt = 1.
+                ref[6] = 1.
+                H = ref[:3]
+                phi,beta = G2lat.CrsAng(H,cell,SGData)
+                psi,gam,x,x = G2lat.SamAng(ref[3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
+                for item in parmDict:
+                    if 'C' in item:
+                        L,M,N = eval(item.strip('C'))
+                        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
+                        Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
+                        Lnorm = G2lat.Lnorm(L)
+                        ref[6] += parmDict[item]*Lnorm*Kcl*Ksl
+                mat.append(wt*(ref[5]-ref[6]))
+            Mat = np.concatenate((Mat,np.array(mat)))
+        return Mat
+        
+    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList):
+        Mat = np.empty(0)
+        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
+        for hist in Gangls.keys():
+            mat = np.zeros((len(refData[hist]),len(varyList)))
+            for i,ref in enumerate(refData[hist]):
+#                wt = np.sqrt(ref[4])
+                wt = 1.
+                H = ref[:3]
+                phi,beta = G2lat.CrsAng(H,cell,SGData)
+                psi,gam,dPdA,dGdA = G2lat.SamAng(ref[3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
+                for j,item in enumerate(varyList):
+                    if 'C' in item:
+                        L,M,N = eval(item.strip('C'))
+                        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
+                        Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
+                        Lnorm = G2lat.Lnorm(L)
+                        mat[i,j] = -wt*Lnorm*Kcl*Ksl
+                        for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
+                            try:
+                                l = varyList.index(itema)
+                                mat[i,l] += parmDict[item]*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
+                            except ValueError:
+                                pass
+            if len(Mat):
+                Mat = np.concatenate((Mat,mat))
+            else:
+                Mat = mat
+        return Mat
+
+    print ' Fit texture for '+General['Name']
     SGData = General['SGData']
+    cell = General['Cell'][1:7]
     Texture = General['SH Texture']
     if not Texture['Order']:
         return 'No spherical harmonics coefficients'
-        
+    varyList = []
+    parmDict = copy.copy(Texture['SH Coeff'][1])
+    for item in ['Sample omega','Sample chi','Sample phi']:
+        parmDict[item] = Texture[item][1]
+        if Texture[item][0]:
+            varyList.append(item)
+    if Texture['SH Coeff'][0]:
+        varyList += Texture['SH Coeff'][1].keys()
+    while True:
+        begin = time.time()
+        values =  np.array(Dict2Values(parmDict, varyList))
+        result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,
+            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList))
+        ncyc = int(result[2]['nfev']/2)
+        if ncyc:
+            runtime = time.time()-begin    
+            chisq = np.sum(result[2]['fvec']**2)
+            Values2Dict(parmDict, varyList, result[0])
+            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
+            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(result[2]['fvec']),' Number of parameters: ',len(varyList)
+            print 'refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
+            print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)
+            try:
+                sig = np.sqrt(np.diag(result[1])*GOF)
+                if np.any(np.isnan(sig)):
+                    print '*** Least squares aborted - some invalid esds possible ***'
+                break                   #refinement succeeded - finish up!
+            except ValueError:          #result[1] is None on singular matrix
+                print '**** Refinement failed - singular matrix ****'
+                return None
+        else:
+            break
+    
+    for hist in refData:
+        print ' Texture corrections for '+hist
+        for ref in refData[hist]:
+            print ' %d %d %d %.3f %.3f'%(int(ref[0]),int(ref[1]),int(ref[2]),ref[5],ref[6])
+    if ncyc:
+        for parm in parmDict:
+            if 'C' in parm:
+                Texture['SH Coeff'][1][parm] = parmDict[parm]
+            else:
+                Texture[parm][1] = parmDict[parm]  
+        sigDict = dict(zip(varyList,sig))
+        printSpHarm(Texture,sigDict)
         
     return None

trunk/GSASIIobj.py ¶

@@ -848 +848 @@
 
 DefaultControls = {
-    'deriv type':'analytic Hessian',    #default controls
-    'min dM/M':0.0001,'shift factor':1.,'max cyc':3,'F**2':True,
-    'minF/sig':0,
+    'deriv type':'analytic Hessian',
+    'min dM/M':0.0001,'shift factor':1.,'max cyc':3,'F**2':False,
+    'UsrReject':{'minF/sig':0,'MinExt':0.01,'MaxDF/F':100.,'MaxD':500.,'MinD':0.05},
+    'Copy2Next':False,'Reverse Seq':False,
     'Author':'no name',
     'FreePrm1':'Sample humidity (%)',
     'FreePrm2':'Sample voltage (V)',
     'FreePrm3':'Applied load (MN)',
+    'SeqPseudoVars':{},'SeqParFitEqList':[],'ShowCell':False,
     }
 '''Values to be used as defaults for the initial contents of the ``Controls``

trunk/GSASIIphsGUI.py ¶

@@ -3516 +3516 @@
 ################################################################################
         
-    def UpdateTexture():        
+    def UpdateTexture():
+                
         def SetSHCoef():
             cofNames = G2lat.GenSHCoeff(SGData['SGLaue'],SamSym[textureData['Model']],textureData['Order'])
@@ -5927 +5928 @@
         histNames = []
         refData = {}
-        SamAngs = {}
+        Gangls = {}
         for name in Histograms.keys():
             if 'PWDR' in name:
@@ -5934 +5935 @@
                 histNames.append(name)
                 Id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,name)
+                Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'))
                 Sample = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Sample Parameters'))
-                SamAngs[name] = copy.copy([Sample[item] for item in['Omega','Chi','Phi','Azimuth']])
+                Gangls[name] = copy.copy([Sample[item] for item in['Phi','Chi','Omega','Azimuth']])
                 RefDict = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Reflection Lists'))[phaseName]
                 Refs = RefDict['RefList'].T  #np.array!
-                if 'T' in RefDict['Type']: it = 3  #TOF offset for alp, bet, wave
                 if RefDict['Super']: im = 1     #(3+1) offset for m
-                refData[name] = np.column_stack((Refs[0],Refs[1],Refs[2],Refs[8+im],Refs[12+im+it]))
-        Error = G2mth.FitTexture(General,SamAngs,refData)
+                if 'T' in RefDict['Type']: 
+                    it = 3  #TOF offset for alp, bet, wave
+                    tth = np.ones_like(Refs[0])*Inst[0]['2-theta'][0]
+                    refData[name] = np.column_stack((Refs[0],Refs[1],Refs[2],tth,Refs[8+im],Refs[12+im+it],np.zeros_like(Refs[0])))
+                else:   # xray - typical caked 2D image data
+                    refData[name] = np.column_stack((Refs[0],Refs[1],Refs[2],Refs[5+im],Refs[8+im],Refs[12+im+it],np.zeros_like(Refs[0])))
+        Error = G2mth.FitTexture(General,Gangls,refData)
         if Error:
             wx.MessageBox(Error,caption='Fit Texture Error',style=wx.ICON_EXCLAMATION)
-        G2ddG.UpdateDData(G2frame,DData,data)
+        UpdateTexture()
+        G2plt.PlotTexture(G2frame,data,Start=False)            
             
     def OnTextureClear(event):

trunk/GSASIIstrMath.py ¶

@@ -1232 +1232 @@
     else:
         tth = refl[5+im]
-    FORPI = 4.0*np.pi
     IFCoup = 'Bragg' in calcControls[hfx+'instType']
     odfCor = 1.0
@@ -1286 +1285 @@
     else:
         tth = refl[5+im]
-    FORPI = 12.5663706143592
     odfCor = 1.0
     dFdODF = {}