Changeset 1790


Ignore:
Timestamp:
Apr 15, 2015 4:20:11 PM (7 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:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1786 r1790  
    3333import GSASIIpwd as G2pwd
    3434import numpy.fft as fft
     35import scipy.optimize as so
    3536import pypowder as pyd
    3637
     
    15921593################################################################################
    15931594
    1594 def FitTexture(General,SamAngs,refData):
    1595     print 'fit texture for '+General["Name"]
     1595def FitTexture(General,Gangls,refData):
     1596   
     1597    def printSpHarm(textureData,SHtextureSig):
     1598        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
     1599        names = ['omega','chi','phi']
     1600        namstr = '  names :'
     1601        ptstr =  '  values:'
     1602        sigstr = '  esds  :'
     1603        for name in names:
     1604            namstr += '%12s'%('Sample '+name)
     1605            ptstr += '%12.3f'%(textureData['Sample '+name][1])
     1606            if 'Sample '+name in SHtextureSig:
     1607                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
     1608            else:
     1609                sigstr += 12*' '
     1610        print namstr
     1611        print ptstr
     1612        print sigstr
     1613        print '\n Texture coefficients:'
     1614        SHcoeff = textureData['SH Coeff'][1]
     1615        SHkeys = SHcoeff.keys()
     1616        nCoeff = len(SHcoeff)
     1617        nBlock = nCoeff/10+1
     1618        iBeg = 0
     1619        iFin = min(iBeg+10,nCoeff)
     1620        for block in range(nBlock):
     1621            namstr = '  names :'
     1622            ptstr =  '  values:'
     1623            sigstr = '  esds  :'
     1624            for name in SHkeys[iBeg:iFin]:
     1625                if 'C' in name:
     1626                    namstr += '%12s'%(name)
     1627                    ptstr += '%12.3f'%(SHcoeff[name])
     1628                    if name in SHtextureSig:
     1629                        sigstr += '%12.3f'%(SHtextureSig[name])
     1630                    else:
     1631                        sigstr += 12*' '
     1632            print namstr
     1633            print ptstr
     1634            print sigstr
     1635            iBeg += 10
     1636            iFin = min(iBeg+10,nCoeff)
     1637           
     1638    def Dict2Values(parmdict, varylist):
     1639        '''Use before call to leastsq to setup list of values for the parameters
     1640        in parmdict, as selected by key in varylist'''
     1641        return [parmdict[key] for key in varylist]
     1642       
     1643    def Values2Dict(parmdict, varylist, values):
     1644        ''' Use after call to leastsq to update the parameter dictionary with
     1645        values corresponding to keys in varylist'''
     1646        parmdict.update(zip(varylist,values))
     1647       
     1648    def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList):
     1649        parmDict.update(zip(varyList,values))
     1650        Mat = np.empty(0)
     1651        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
     1652        for hist in Gangls.keys():
     1653            mat = []
     1654            for ref in refData[hist]:
     1655#                wt = np.sqrt(ref[4])
     1656                wt = 1.
     1657                ref[6] = 1.
     1658                H = ref[:3]
     1659                phi,beta = G2lat.CrsAng(H,cell,SGData)
     1660                psi,gam,x,x = G2lat.SamAng(ref[3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
     1661                for item in parmDict:
     1662                    if 'C' in item:
     1663                        L,M,N = eval(item.strip('C'))
     1664                        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
     1665                        Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
     1666                        Lnorm = G2lat.Lnorm(L)
     1667                        ref[6] += parmDict[item]*Lnorm*Kcl*Ksl
     1668                mat.append(wt*(ref[5]-ref[6]))
     1669            Mat = np.concatenate((Mat,np.array(mat)))
     1670        return Mat
     1671       
     1672    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList):
     1673        Mat = np.empty(0)
     1674        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
     1675        for hist in Gangls.keys():
     1676            mat = np.zeros((len(refData[hist]),len(varyList)))
     1677            for i,ref in enumerate(refData[hist]):
     1678#                wt = np.sqrt(ref[4])
     1679                wt = 1.
     1680                H = ref[:3]
     1681                phi,beta = G2lat.CrsAng(H,cell,SGData)
     1682                psi,gam,dPdA,dGdA = G2lat.SamAng(ref[3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
     1683                for j,item in enumerate(varyList):
     1684                    if 'C' in item:
     1685                        L,M,N = eval(item.strip('C'))
     1686                        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
     1687                        Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
     1688                        Lnorm = G2lat.Lnorm(L)
     1689                        mat[i,j] = -wt*Lnorm*Kcl*Ksl
     1690                        for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
     1691                            try:
     1692                                l = varyList.index(itema)
     1693                                mat[i,l] += parmDict[item]*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
     1694                            except ValueError:
     1695                                pass
     1696            if len(Mat):
     1697                Mat = np.concatenate((Mat,mat))
     1698            else:
     1699                Mat = mat
     1700        return Mat
     1701
     1702    print ' Fit texture for '+General['Name']
    15961703    SGData = General['SGData']
     1704    cell = General['Cell'][1:7]
    15971705    Texture = General['SH Texture']
    15981706    if not Texture['Order']:
    15991707        return 'No spherical harmonics coefficients'
    1600        
     1708    varyList = []
     1709    parmDict = copy.copy(Texture['SH Coeff'][1])
     1710    for item in ['Sample omega','Sample chi','Sample phi']:
     1711        parmDict[item] = Texture[item][1]
     1712        if Texture[item][0]:
     1713            varyList.append(item)
     1714    if Texture['SH Coeff'][0]:
     1715        varyList += Texture['SH Coeff'][1].keys()
     1716    while True:
     1717        begin = time.time()
     1718        values =  np.array(Dict2Values(parmDict, varyList))
     1719        result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,
     1720            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList))
     1721        ncyc = int(result[2]['nfev']/2)
     1722        if ncyc:
     1723            runtime = time.time()-begin   
     1724            chisq = np.sum(result[2]['fvec']**2)
     1725            Values2Dict(parmDict, varyList, result[0])
     1726            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
     1727            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(result[2]['fvec']),' Number of parameters: ',len(varyList)
     1728            print 'refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
     1729            print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)
     1730            try:
     1731                sig = np.sqrt(np.diag(result[1])*GOF)
     1732                if np.any(np.isnan(sig)):
     1733                    print '*** Least squares aborted - some invalid esds possible ***'
     1734                break                   #refinement succeeded - finish up!
     1735            except ValueError:          #result[1] is None on singular matrix
     1736                print '**** Refinement failed - singular matrix ****'
     1737                return None
     1738        else:
     1739            break
     1740   
     1741    for hist in refData:
     1742        print ' Texture corrections for '+hist
     1743        for ref in refData[hist]:
     1744            print ' %d %d %d %.3f %.3f'%(int(ref[0]),int(ref[1]),int(ref[2]),ref[5],ref[6])
     1745    if ncyc:
     1746        for parm in parmDict:
     1747            if 'C' in parm:
     1748                Texture['SH Coeff'][1][parm] = parmDict[parm]
     1749            else:
     1750                Texture[parm][1] = parmDict[parm] 
     1751        sigDict = dict(zip(varyList,sig))
     1752        printSpHarm(Texture,sigDict)
    16011753       
    16021754    return None
  • trunk/GSASIIobj.py

    r1698 r1790  
    848848
    849849DefaultControls = {
    850     'deriv type':'analytic Hessian',    #default controls
    851     'min dM/M':0.0001,'shift factor':1.,'max cyc':3,'F**2':True,
    852     'minF/sig':0,
     850    'deriv type':'analytic Hessian',
     851    'min dM/M':0.0001,'shift factor':1.,'max cyc':3,'F**2':False,
     852    'UsrReject':{'minF/sig':0,'MinExt':0.01,'MaxDF/F':100.,'MaxD':500.,'MinD':0.05},
     853    'Copy2Next':False,'Reverse Seq':False,
    853854    'Author':'no name',
    854855    'FreePrm1':'Sample humidity (%)',
    855856    'FreePrm2':'Sample voltage (V)',
    856857    'FreePrm3':'Applied load (MN)',
     858    'SeqPseudoVars':{},'SeqParFitEqList':[],'ShowCell':False,
    857859    }
    858860'''Values to be used as defaults for the initial contents of the ``Controls``
  • trunk/GSASIIphsGUI.py

    r1785 r1790  
    35163516################################################################################
    35173517       
    3518     def UpdateTexture():       
     3518    def UpdateTexture():
     3519               
    35193520        def SetSHCoef():
    35203521            cofNames = G2lat.GenSHCoeff(SGData['SGLaue'],SamSym[textureData['Model']],textureData['Order'])
     
    59275928        histNames = []
    59285929        refData = {}
    5929         SamAngs = {}
     5930        Gangls = {}
    59305931        for name in Histograms.keys():
    59315932            if 'PWDR' in name:
     
    59345935                histNames.append(name)
    59355936                Id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,name)
     5937                Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'))
    59365938                Sample = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Sample Parameters'))
    5937                 SamAngs[name] = copy.copy([Sample[item] for item in['Omega','Chi','Phi','Azimuth']])
     5939                Gangls[name] = copy.copy([Sample[item] for item in['Phi','Chi','Omega','Azimuth']])
    59385940                RefDict = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Reflection Lists'))[phaseName]
    59395941                Refs = RefDict['RefList'].T  #np.array!
    5940                 if 'T' in RefDict['Type']: it = 3  #TOF offset for alp, bet, wave
    59415942                if RefDict['Super']: im = 1     #(3+1) offset for m
    5942                 refData[name] = np.column_stack((Refs[0],Refs[1],Refs[2],Refs[8+im],Refs[12+im+it]))
    5943         Error = G2mth.FitTexture(General,SamAngs,refData)
     5943                if 'T' in RefDict['Type']:
     5944                    it = 3  #TOF offset for alp, bet, wave
     5945                    tth = np.ones_like(Refs[0])*Inst[0]['2-theta'][0]
     5946                    refData[name] = np.column_stack((Refs[0],Refs[1],Refs[2],tth,Refs[8+im],Refs[12+im+it],np.zeros_like(Refs[0])))
     5947                else:   # xray - typical caked 2D image data
     5948                    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])))
     5949        Error = G2mth.FitTexture(General,Gangls,refData)
    59445950        if Error:
    59455951            wx.MessageBox(Error,caption='Fit Texture Error',style=wx.ICON_EXCLAMATION)
    5946         G2ddG.UpdateDData(G2frame,DData,data)
     5952        UpdateTexture()
     5953        G2plt.PlotTexture(G2frame,data,Start=False)           
    59475954           
    59485955    def OnTextureClear(event):
  • trunk/GSASIIstrMath.py

    r1787 r1790  
    12321232    else:
    12331233        tth = refl[5+im]
    1234     FORPI = 4.0*np.pi
    12351234    IFCoup = 'Bragg' in calcControls[hfx+'instType']
    12361235    odfCor = 1.0
     
    12861285    else:
    12871286        tth = refl[5+im]
    1288     FORPI = 12.5663706143592
    12891287    odfCor = 1.0
    12901288    dFdODF = {}
Note: See TracChangeset for help on using the changeset viewer.