Changeset 657


Ignore:
Timestamp:
Jun 26, 2012 3:44:56 PM (9 years ago)
Author:
vondreele
Message:

fixups of findOffset & mapsearch

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r655 r657  
    610610## keep this
    611611               
    612 def findOffset(SGData,rho,Fhkl):
    613    
    614     def calcPhase(DX,DH,Dphi):
    615         H,K,L = DH.T
    616         Phi = (H*DX[0]+K*DX[1]+L*DX[2]+0.5) % 1.
    617         return Dphi-Phi
    618    
     612def findOffset(SGData,rho,Fhkl):   
    619613    if SGData['SpGrp'] == 'P 1':
    620614        return [0,0,0]   
    621 # will need to consider 'SGPolax': one of '','x','y','x y','z','x z','y z','xyz','111'
    622615    mapShape = rho.shape
    623616    steps = np.array(mapShape)
     
    660653    DH = np.array(DH)
    661654    Dphi = np.array(Dphi)
    662 #    for item in zip(DH,Dphi):
    663 #        print item[0],'%.4f'%(item[1])
    664     DX = np.zeros(3)
    665     X,Y,Z = np.mgrid[0:1:10j,0:1:10j,0:1:10j]
     655    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
    666656    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
    667     Mmin = 1.e10
    668    
    669     for xyz in XYZ:             #do a global search for best roll
    670         M = np.sum(calcPhase(xyz,DH,Dphi)**2)
    671         if M < Mmin:
    672             DX = xyz
    673             Mmin = M
    674    
    675     result = so.leastsq(calcPhase,DX,full_output=True,args=(DH,Dphi))
    676     for item in zip(DH,Dphi,result[2]['fvec']):
    677         print item[0],'%.4f %.4f'%(item[1],item[2])
    678     chisq = np.sum(result[2]['fvec']**2)
    679     DX = np.array(np.fix(-result[0]*steps),dtype='i')
     657    Mmap = np.reshape(np.sum(((np.dot(XYZ,DH.T)+.5)%1.-Dphi)**2,axis=1),newshape=steps)
     658    chisq = np.min(Mmap)
     659    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
    680660    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
    681661    return DX
    682662   
    683 def findOffset2(SGData,rho,Fhkl):
    684     if SGData['SpGrp'] == 'P 1':
    685         return [0,0,0]
    686     mapShape = rho.shape
    687     mapHalf = np.array(mapShape)/2
    688     steps = np.array(mapShape)
    689     hklShape = Fhkl.shape
    690     hklHalf = np.array(hklShape)/2
    691     for M,T in SGData['SGOps'][1:]:
    692         FThkl = np.zeros(shape=hklShape,dtype='c16')
    693         for ind,F in enumerate(Fhkl.flatten()):
    694             if F:
    695                 HKL = np.unravel_index(ind,hklShape)-hklHalf
    696                 HKLT = np.inner(HKL,M.T)
    697                 phi = (np.angle(F,deg=True)/360.+np.vdot(T,HKL) % 1.)*360.
    698                 phasep = complex(cosd(phi),-sind(phi))      #want F*
    699                 h,k,l = HKLT+hklHalf
    700                 FThkl[h,k,l] = np.absolute(F)*phasep
    701         Cd = np.real(fft.fftn(fft.fftshift(Fhkl*FThkl)))
    702         CdMax = np.array(np.unravel_index(np.argmax(Cd),mapShape))
    703         print CdMax,CdMax/2
    704     DX = CdMax/2           
    705    
    706     print ' Test map offset : %d %d %d'%(DX[0],DX[1],DX[2])
    707     return DX
    708 
    709 
    710663def ChargeFlip(data,reflData,pgbar):
    711664    generalData = data['General']
     
    773726                break
    774727            GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
    775             if not GoOn:
     728            if not GoOn or Ncyc > 10000:
    776729                break
    777730        except FloatingPointError:
     
    782735    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
    783736    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
    784 #    roll = findOffset2(SGData,CErho,CEhkl) #doesn't work
    785737    roll = findOffset(SGData,CErho,CEhkl)
    786738       
     
    791743    return mapData
    792744   
    793 def SearchMap(data,keepDup=False):
     745def SearchMap(data,keepDup=False,Pgbar=None):
    794746   
    795747    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
     
    811763    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
    812764        Mag,x0,y0,z0,sig = parms
    813         if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
    814             return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(x0-rX)*(y0-rY)+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
    815         else:
    816             return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
     765#        if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
     766#            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(x0-rX)*(y0-rY)+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
     767#        else:
     768#            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
     769        return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
    817770       
    818771    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
     
    887840                peaks.append(peak)
    888841                mags.append(x1[0])
    889             if len(peaks) > 300:
     842            GoOn = Pgbar.Update(len(peaks),newmsg='%s%d'%('No. Peaks found =',len(peaks)))[0]
     843            if not GoOn or len(peaks) > 300:
    890844                break
    891845        rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(x1,rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
  • trunk/GSASIIphsGUI.py

    r652 r657  
    38253825        print ' Begin fourier map search - can take some time'
    38263826        time0 = time.time()
    3827         wx.BeginBusyCursor()
     3827
     3828        pgbar = wx.ProgressDialog('Map search','No. Peaks found =',301.0,
     3829            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_CAN_ABORT)
     3830        screenSize = wx.ClientDisplayRect()
     3831        Size = pgbar.GetSize()
     3832        Size = (int(Size[0]*1.2),Size[1]) # increase size a bit along x
     3833        pgbar.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
     3834        pgbar.SetSize(Size)
    38283835        try:
    3829             peaks,mags = G2mth.SearchMap(data,keepDup=True)
     3836            peaks,mags = G2mth.SearchMap(data,keepDup=True,Pgbar=pgbar)
    38303837        finally:
    3831             wx.EndBusyCursor()
     3838            pgbar.Destroy()
     3839
     3840
     3841#        wx.BeginBusyCursor()
     3842#        try:
     3843#            peaks,mags = G2mth.SearchMap(data,keepDup=True)
     3844#        finally:
     3845#            wx.EndBusyCursor()
    38323846        sortIdx = np.argsort(mags.flatten())
    38333847        if len(peaks):
  • trunk/GSASIIstruct.py

    r656 r657  
    525525                if isum == jsum:
    526526                    eqvDict[varyI].append(varyJ)
    527             elif abs(dspI-dspJ) < 1.e-4:
     527            elif abs(dspI-dspJ)/dspI < 1.e-4:
    528528                eqvDict[varyI].append(varyJ)
    529529    for item in pawleyVary:
Note: See TracChangeset for help on using the changeset viewer.