Changeset 648


Ignore:
Timestamp:
Jun 15, 2012 1:57:29 AM (9 years ago)
Author:
vondreele
Message:

1st stab at a unique atoms picker - mostly GUI
format fix for constraint variable lister
fix to findOffset to off by one problem
make multiplicities in PWDR & HKLF sets jive
put Lorentz factor in Pawley estimator & add FWHM calculator
fix to instprm output
more work on Pawley penalty fxn - closer!!

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIgrid.py

    r635 r648  
    3737htmlFirstUse = True
    3838
    39 [ wxID_FOURCALC,wxID_FOURSEARCH, wxID_PEAKSMOVE, wxID_PEAKSCLEAR, wxID_CHARGEFLIP,
    40 ] = [wx.NewId() for item in range(5)]
     39[ wxID_FOURCALC, wxID_FOURSEARCH, wxID_PEAKSMOVE, wxID_PEAKSCLEAR, wxID_CHARGEFLIP,
     40    wxID_PEAKSUNIQUE,
     41] = [wx.NewId() for item in range(6)]
    4142
    4243[ wxID_PWDRADD, wxID_HKLFADD, wxID_DATADELETE,
     
    565566        self.MapPeaksEdit.Append(id=wxID_PEAKSMOVE, kind=wx.ITEM_NORMAL,text='Move peaks',
    566567            help='Move selected peaks to atom list')
     568        self.MapPeaksEdit.Append(id=wxID_PEAKSUNIQUE, kind=wx.ITEM_NORMAL,text='Unique peaks',
     569            help='Reduce map peak list to unique set')
    567570        self.MapPeaksEdit.Append(id=wxID_PEAKSCLEAR, kind=wx.ITEM_NORMAL,text='Clear peaks',
    568571            help='Clear the map peak list')
  • trunk/GSASIImapvars.py

    r584 r648  
    677677    printlist.append(3*[None])
    678678    for name,val,esd in printlist:
    679         if len(s1) > 110 or name is None:
     679        if len(s1) > 120 or name is None:
    680680            print
    681681            print s1
     
    688688            s2 = ' value :'
    689689            s3 = ' sig   :'
    690         s1 += '%12s' % (name)
    691         s2 += '%12.6f' % (val)
     690        s1 += '%15s' % (name)
     691        s2 += '%15.4f' % (val)
    692692        if esd is None:
    693             s3 += '%12s' % ('n/a')
     693            s3 += '%15s' % ('n/a')
    694694        else:   
    695             s3 += '%12.6f' % (esd)
     695            s3 += '%15.4f' % (esd)
    696696
    697697def ComputeDepESD(covMatrix,varyList,parmDict):
  • trunk/GSASIImath.py

    r647 r648  
    650650#        print item[0],'%.4f %.4f'%(item[1],item[2])
    651651    chisq = np.sum(result[2]['fvec']**2)
    652     DX = np.array(np.rint(-result[0]*steps),dtype='i')
     652    DX = np.array(np.fix(-result[0]*steps),dtype='i')
    653653    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
    654654    return DX
     
    749749    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
    750750   
    751     def noDuplicate(xyz,peaks,SGData):                  #be careful where this is used - it's slow
     751    def noDuplicate(xyz,peaks,SGData,incre):                  #be careful where this is used - it's slow
    752752        equivs = G2spc.GenAtom(xyz,SGData)
    753753        xyzs = [equiv[0] for equiv in equivs]
    754754        for x in xyzs:
    755             if True in [np.allclose(x,peak,atol=0.002) for peak in peaks]:
     755            if True in [np.allclose(x,peak,atol=0.02) for peak in peaks]:
    756756                return False
    757757        return True
     
    838838                peaks.append(peak)
    839839                mags.append(x1[0])
    840             elif noDuplicate(peak,peaks,SGData) and x1[0] > 0.:
     840            elif noDuplicate(peak,peaks,SGData,incre) and x1[0] > 0.:
    841841                peaks.append(peak)
    842842                mags.append(x1[0])
  • trunk/GSASIIphsGUI.py

    r647 r648  
    36223622                ext,mul = G2spc.GenHKLf([h,k,l],SGData)[:2]
    36233623                if not ext:
     3624                    mul *= 2        #for powder multiplicity
    36243625                    PawleyPeaks.append([h,k,l,mul,d,False,100.0,1.0])
    36253626        finally:
     
    36573658            indx = np.searchsorted(xdata[0],pos)
    36583659            try:
    3659                 ref[6] = xdata[1][indx]/ref[3]
    3660                 pola,dIdPola = G2pwd.Polarization(Inst['Polariz.'],xdata[0][indx],0.0)
    3661                 ref[6] /= pola
     3660                FWHM = max(0.001,G2pwd.getFWHM(pos,Inst))/2.
     3661                dx = xdata[0][indx+1]-xdata[0][indx]
     3662                ref[6] = FWHM*xdata[1][indx]/dx
     3663                Lorenz = 1./(2.*sind(xdata[0][indx]/2.)**2*cosd(xdata[0][indx]/2.))           #Lorentz correction
     3664                pola,dIdPola = G2pwd.Polarization(Inst['Polariz.'],xdata[0][indx],Inst['Azimuth'])
     3665                ref[6] /= (Lorenz*pola*ref[3])
    36623666            except IndexError:
    36633667                pass
     
    37333737        FillMapPeaksGrid()
    37343738        G2plt.PlotStructure(G2frame,data)
     3739       
     3740    def OnPeaksUnique(event):
     3741        generalData = data['General']
     3742        SGData = generalData['SGData']
     3743        if 'Map Peaks' in data:
     3744            mapPeaks = data['Map Peaks']
     3745            for ipk,peak in enumerate(mapPeaks):
     3746                XYZ = peak[1:]                       
     3747                Equiv = G2spc.GenAtom(XYZ,SGData,Move=True)[1:]     #remove self
    37353748   
    37363749    def OnFourierMaps(event):
     
    38083821        G2frame.dataFrame.SetMenuBar(G2frame.dataFrame.MapPeaksMenu)
    38093822        G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksMove, id=G2gd.wxID_PEAKSMOVE)
     3823        G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksUnique, id=G2gd.wxID_PEAKSUNIQUE)
    38103824        G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksClear, id=G2gd.wxID_PEAKSCLEAR)
    38113825        UpdateDrawAtoms()
  • trunk/GSASIIpwd.py

    r482 r648  
    133133            which (if either) of these is "right"?
    134134    return:
    135         pola = (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
    136             Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
     135        pola = ((1-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \
     136            (1-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
    137137        dpdPola: derivative needed for least squares
    138138    """
     
    585585        fmin,fmax = [fmax,fmin]         
    586586    return widths,fmin,fmax
     587   
     588def getFWHM(TTh,Inst):
     589    sig = lambda Th,U,V,W: 1.17741*math.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))*math.pi/180.
     590    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/180.
     591    gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
     592    s = sig(TTh/2.,Inst['U'],Inst['V'],Inst['W'])*100.
     593    g = gam(TTh/2.,Inst['X'],Inst['Y'])*100.
     594    print TTh,s,g
     595    return gamFW(g,s)
     596   
    587597               
    588598def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
  • trunk/GSASIIpwdGUI.py

    r645 r648  
    680680                filename = dlg.GetPath()
    681681                File = open(filename,'w')
    682                 File.write("#GSAS-II instrument parameter file; do not add/delete or change order of items!")
     682                File.write("#GSAS-II instrument parameter file; do not add/delete or change order of items!\n")
    683683                for i,item in enumerate(data[3]):
    684684                    File.write(item+':'+str(data[1][i])+'\n')
  • trunk/GSASIIstruct.py

    r646 r648  
    491491       
    492492def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
    493     if SGLaue in ['-1','2/m','mmm']:
    494         return                      #no Pawley symmetry required constraints
     493#    if SGLaue in ['-1','2/m','mmm']:
     494#        return                      #no Pawley symmetry required constraints
    495495    eqvDict = {}
    496496    for i,varyI in enumerate(pawleyVary):
     
    525525                if isum == jsum:
    526526                    eqvDict[varyI].append(varyJ)
    527 #            if abs(dspI-dspJ) < 3.e-4:
    528 #                eqvDict[varyI].append(varyJ)
     527            elif abs(dspI-dspJ) < 1.e-4:
     528                eqvDict[varyI].append(varyJ)
    529529    for item in pawleyVary:
    530530        if eqvDict[item]:
     
    18361836    for i,item in enumerate(varyList):
    18371837        if 'PWLref' in item and parmDict[item] < 0.:
    1838             pFxn[i] = -parmDict[item]
     1838            pFxn[i] = -parmDict[item]**2        #checked OK
    18391839    return pFxn
    18401840   
     
    18431843    for i,item in enumerate(varyList):
    18441844        if 'PWLref' in item and parmDict[item] < 0.:
    1845             pDerv[i] += -1./parmDict[item]
     1845            pDerv[i] += 2.*parmDict[item]
    18461846    return pDerv
    18471847
     
    24672467            xF = np.searchsorted(x,Limits[1])
    24682468            ymb = np.array(y-yb)
     2469            ymb = np.where(ymb==0.,1.0,ymb)
    24692470            ycmb = np.array(yc-yb)
    2470             ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
     2471            ratio = ymb/ycmb           
    24712472            refLists = Histogram['Reflection Lists']
    24722473            for phase in refLists:
     
    25712572                if Phase['General'].get('doPawley'):
    25722573                    try:
    2573                         refl[9] = parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]
     2574                       refl[9] = parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]
    25742575                    except KeyError:
    25752576#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
     
    27192720                        if Ka2:
    27202721                            dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9]
    2721                         # Assuming Pawley variables not in constraints
    27222722                    except ValueError:
    27232723                        pass
     
    28922892            dMdv = dMdvh
    28932893           
    2894 #    dpdv = penaltyDeriv(parmdict,varylist)
    2895 #    dMdv = np.concatenate((dMdv.T,np.outer(dpdv,dpdv))).T
     2894    dpdv = penaltyDeriv(parmdict,varylist)
     2895    if np.any(dpdv):
     2896        dMdv = np.concatenate((dMdv.T,np.outer(dpdv,dpdv))).T
    28962897    return dMdv
    28972898
     
    29712972        else:
    29722973            continue        #skip non-histogram entries
    2973 
    2974 #    dpdv = penaltyDeriv(parmdict,varylist)
    2975 #    Vec += dpdv*penaltyFxn(parmdict,varylist)
    2976 #    Hess += np.outer(dpdv,dpdv)
     2974    dpdv = penaltyDeriv(parmdict,varylist)
     2975    if np.any(dpdv):
     2976        Vec += dpdv*penaltyFxn(parmdict,varylist)
     2977        pHess = np.zeros((len(varylist),len(varylist)))
     2978        for i,val in enumerate(dpdv):
     2979            pHess[i][i] = dpdv[i]**2
     2980        Hess += pHess
    29772981    return Vec,Hess
    29782982
     
    30793083            parmDict['saved values'] = values
    30803084            raise Exception         #Abort!!
    3081 #    M = np.concatenate((M,penaltyFxn(parmdict,varylist)))
     3085    pFunc = penaltyFxn(parmdict,varylist)
     3086    if np.any(pFunc):
     3087        M = np.concatenate((M,pFunc))
    30823088    return M
    30833089                       
     
    31543160                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
    31553161            ncyc = result[2]['num cyc']+1
    3156             Rvals['lamMax'] = result[2]['lamMax']                           
     3162            Rvals['lamMax'] = result[2]['lamMax']
    31573163        else:           #'numeric'
    31583164            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
Note: See TracChangeset for help on using the changeset viewer.