Changeset 530


Ignore:
Timestamp:
Apr 6, 2012 1:30:45 PM (10 years ago)
Author:
vondreele
Message:

small fix to HessianLSQ about # of cycles
finish map peak search (no GUI/plot output yet)

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r486 r530  
    9191    lamMax = lam
    9292    nfev = 0
    93     while icycle <= maxcyc:
     93    while icycle < maxcyc:
    9494        lamMax = max(lamMax,lam)
    9595        M = func(x0,*args)
     
    131131        Bmat = nl.inv(Amat)
    132132        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
    133     except LinAlgError:
    134         psing = list(np.where(np.diag(nl.gr(Amat)[1]) < 1.e-14)[0])
     133    except nl.LinAlgError:
     134        psing = []
     135        if maxcyc:
     136            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
    135137        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
    136 
    137 def calcFouriermap():
    138     print 'Calculate Fourier map'
    139 
    140    
     138   
    141139def getVCov(varyNames,varyList,covMatrix):
    142140    vcov = np.zeros((len(varyNames),len(varyNames)))
  • trunk/GSASIIphsGUI.py

    r529 r530  
    2525import GSASIIIO as G2IO
    2626import GSASIIstruct as G2str
     27import GSASIImath as G2mth
    2728import numpy as np
    2829import numpy.linalg as nl
     
    34453446        data['Drawing']['mapSize'] = 10.
    34463447        print mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f'%(np.max(mapData['rho']),np.min(mapData['rho']))
    3447         G2plt.PlotStructure(G2frame,data)                   
    3448 ## map printing for testing purposes
    3449 #        ix,jy,kz = mapData['rho'].shape
    3450 #        for k in range(kz):
    3451 #            print 'k = ',k
    3452 #            for j in range(jy):
    3453 #                line = ''
    3454 #                if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
    3455 #                    line += (jy-j)*'  '
    3456 #                for i in range(ix):
    3457 #                    r = int(100*mapData['rho'][i,j,k]/mapData['rhoMax'])
    3458 #                    line += '%4d'%(r)
    3459 #                print line+'\n'
    3460 ### keep this               
     3448        G2plt.PlotStructure(G2frame,data)
     3449       
     3450    def printRho(SGData,rho,rhoMax):                         
     3451# map printing for testing purposes
     3452        ix,jy,kz = rho.shape
     3453        for k in range(kz):
     3454            print 'k = ',k
     3455            for j in range(jy):
     3456                line = ''
     3457                if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
     3458                    line += (jy-j)*'  '
     3459                for i in range(ix):
     3460                    r = int(100*rho[i,j,k]/rhoMax)
     3461                    line += '%4d'%(r)
     3462                print line+'\n'
     3463## keep this               
    34613464   
    34623465    def OnSearchMaps(event):
    34633466       
     3467        norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
     3468       
     3469        def noDuplicate(xyz,peaks,SGData):                  #be careful where this is used - it's slow
     3470            equivs = G2spc.GenAtom(xyz,SGData)
     3471            xyzs = [equiv[0] for equiv in equivs]
     3472            for x in xyzs:
     3473                if True in [np.allclose(x,peak,atol=0.002) for peak in peaks]:
     3474                    return False
     3475            return True
     3476               
    34643477        def findRoll(rhoMask,mapHalf):
    34653478            indx = np.array(ma.nonzero(rhoMask)).T
     
    34673480            rhoMax = np.max(rhoList)
    34683481            return mapHalf-indx[np.argmax(rhoList)]
     3482           
     3483        def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
     3484            Mag,x0,y0,z0,sig = parms
     3485            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
     3486                return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(x0-rX)*(y0-rY)+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
     3487            else:
     3488                return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
     3489           
     3490        def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
     3491            Mag,x0,y0,z0,sig = parms
     3492            M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     3493            return M
     3494           
     3495        def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
     3496            Mag,x0,y0,z0,sig = parms
     3497            dMdv = np.zeros(([5,]+list(rX.shape)))
     3498            delt = .01
     3499            for i in range(5):
     3500                parms[i] -= delt
     3501                rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     3502                parms[i] += 2.*delt
     3503                rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     3504                parms[i] -= delt
     3505                dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
     3506            rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     3507            Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
     3508            dMdv = np.reshape(dMdv,(5,rX.size))
     3509            Hess = np.inner(dMdv,dMdv)
     3510           
     3511            return Vec,Hess
    34693512           
    34703513        generalData = data['General']
     
    34793522            rho = copy.copy(mapData['rho'])     #don't mess up original
    34803523            mapHalf = np.array(rho.shape)/2
    3481             step = max(1.0,1./mapData['Resolution'])
    3482             steps = 3*np.array(3*[step,])
     3524            res = mapData['Resolution']
     3525            incre = 1./np.array(rho.shape)
     3526            step = max(1.0,1./res)+1
     3527            steps = np.array(3*[step,])
    34833528        except KeyError:
    34843529            print '**** ERROR - Fourier map not defined'
    34853530            return
    34863531        peaks = []
    3487         while True:
    3488             rhoMask = ma.array(rho,mask=(rho<contLevel))
    3489             rMI = findRoll(rhoMask,mapHalf)
    3490             rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2)
    3491             rMM = mapHalf-steps
    3492             rMP = mapHalf+steps
    3493            
    3494 
    3495             rho = np.roll(np.roll(np.roll(rho,-rMI[2],axis=2),-rMI[1],axis=1),-rMI[0],axis=0)
    3496             break
    3497         print 'search Fourier map'
     3532        mags = []
     3533        print ' Begin fourier map search - can take some time'
     3534        time0 = time.time()
     3535        wx.BeginBusyCursor()
     3536        try:
     3537            while True:
     3538                rhoMask = ma.array(rho,mask=(rho<contLevel))
     3539                if not ma.count(rhoMask):
     3540                    break
     3541                rMI = findRoll(rhoMask,mapHalf)
     3542                rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2)
     3543                rMM = mapHalf-steps
     3544                rMP = mapHalf+steps+1
     3545                rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
     3546                peakInt = np.sum(rhoPeak)*res**3
     3547                rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
     3548                x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
     3549                result = G2mth.HessianLSQ(peakFunc,x0,Hess=peakHess,
     3550                    args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.0001,maxcyc=100)
     3551                x1 = result[0]
     3552                if np.any(x1 < 0):
     3553                    break
     3554                peak = (np.array(x1[1:4])-rMI)*incre
     3555                if not len(peaks):
     3556                    peaks.append(peak)
     3557                    mags.append(x1[0])
     3558                else:
     3559                    if noDuplicate(peak,peaks,SGData) and x1[0] > 0.:
     3560                        peaks.append(peak)
     3561                        mags.append(x1[0])
     3562                rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(result[0],rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
     3563                rho = np.roll(np.roll(np.roll(rho,-rMI[2],axis=2),-rMI[1],axis=1),-rMI[0],axis=0)
     3564        finally:
     3565            wx.EndBusyCursor()
     3566        sortIdx = np.argsort(mags)       
     3567        print ' Map search peaks found:'
     3568        print '  No.    Mag.      x        y        z'
     3569        for j,i in enumerate(sortIdx[::-1]):
     3570            print ' %3d %8.3f %8.5f %8.5f %8.5f'%(j,mags[i],peaks[i][0],peaks[i][1],peaks[i][2])
     3571        print ' Map search finished, time = %.2fs'%(time.time()-time0)
    34983572       
    34993573    def OnTextureRefine(event):
  • trunk/GSASIIspc.py

    r504 r530  
    263263        All  = True return all equivalent positions including duplicates
    264264             = False return only unique positions
    265         Uij  = [U11,U22,U33,U12,U13,U23] or [] if no Uij       
     265        Uij  = [U11,U22,U33,U12,U13,U23] or [] if no Uij
     266        Move = True move generated atom positions to be inside cell
     267             = False do not move atoms       
    266268    return: [[XYZEquiv],Idup,[UijEquiv]]
    267269        [XYZEquiv] is list of equivalent positions (XYZ is first entry)
Note: See TracChangeset for help on using the changeset viewer.