Changeset 949


Ignore:
Timestamp:
Jun 12, 2013 3:08:54 PM (8 years ago)
Author:
vondreele
Message:

changes to MC/SA stuff

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIlattice.py

    r939 r949  
    9393
    9494    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
    95     :param bool inverse: if True return bot G and g; else just G
     95    :param bool inverse: if True return both G and g; else just G
    9696    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
    9797
     
    271271   
    272272def CosSinAngle(U,V,G):
    273     """ calculate sin & cos of angle betwee U & V in generalized coordinates
     273    """ calculate sin & cos of angle between U & V in generalized coordinates
    274274    defined by metric tensor G
    275275
  • trunk/GSASIImath.py

    r948 r949  
    14651465    return Ind
    14661466   
     1467################################################################################
     1468##### single peak fitting profile fxn stuff
     1469################################################################################
     1470
    14671471def getCWsig(ins,pos):
    14681472    tp = tand(pos/2.0)
     
    16281632
    16291633    def update_guess(self, x0):
    1630         pass
     1634        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
    16311635
    16321636    def update_temp(self, x0):
     
    16441648        self.c = self.m * exp(-self.n * self.quench)
    16451649
    1646     def update_guess(self, x0):
    1647         x0 = asarray(x0)
    1648         u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
    1649         T = self.T
    1650         y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
    1651         xc = y*(self.upper - self.lower)
    1652         xnew = x0 + xc
    1653         return xnew
     1650#    def update_guess(self, x0):
     1651#        x0 = asarray(x0)
     1652#        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
     1653#        T = self.T
     1654#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
     1655#        xc = y*(self.upper - self.lower)
     1656#        xnew = x0 + xc
     1657#        return xnew
    16541658
    16551659    def update_temp(self):
     
    16591663
    16601664class cauchy_sa(base_schedule):
    1661     def update_guess(self, x0):
    1662         x0 = asarray(x0)
    1663         numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
    1664         xc = self.learn_rate * self.T * tan(numbers)
    1665         xnew = x0 + xc
    1666         return xnew
     1665#    def update_guess(self, x0):
     1666#        x0 = asarray(x0)
     1667#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
     1668#        xc = self.learn_rate * self.T * tan(numbers)
     1669#        xnew = x0 + xc
     1670#        return xnew
    16671671
    16681672    def update_temp(self):
     
    16721676
    16731677class boltzmann_sa(base_schedule):
    1674     def update_guess(self, x0):
    1675         std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
    1676         x0 = asarray(x0)
    1677         xc = squeeze(random.normal(0, 1.0, size=self.dims))
    1678 
    1679         xnew = x0 + xc*std*self.learn_rate
    1680         return xnew
     1678#    def update_guess(self, x0):
     1679#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
     1680#        x0 = asarray(x0)
     1681#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
     1682#
     1683#        xnew = x0 + xc*std*self.learn_rate
     1684#        return xnew
    16811685
    16821686    def update_temp(self):
     
    16901694        self.__dict__.update(options)
    16911695       
    1692     def update_guess(self,x0):
    1693         x0 = np.asarray(x0)
    1694         u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
    1695         xnew = x0+u
    1696         return xnew
     1696#    def update_guess(self,x0):
     1697#        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
    16971698       
    16981699    def update_temp(self):
    16991700        self.k += 1
    1700         self.T = self.T0*self.slope**k
     1701        self.T = self.T0*self.slope**self.k
    17011702       
    17021703class Tremayne_sa(base_schedule):   #needs fixing for two steps
     
    17051706        self.__dict__.update(options)
    17061707
    1707     def update_guess(self,x0):
    1708         x0 = np.asarray(x0)
    1709         u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
    1710         xnew = x0+u
    1711         return xnew       
     1708#    def update_guess(self,x0):
     1709#        x0 = np.asarray(x0)
     1710#        u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
     1711#        xnew = x0+u
     1712#        return xnew       
    17121713   
    17131714    def update_temp(self):
     
    17301731           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
    17311732           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
    1732            lower=-100, upper=100, dwell=50, slope=0.9):
     1733           lower=-100, upper=100, dwell=50, slope=0.9,dlg=None):
    17331734    """Minimize a function using simulated annealing.
    17341735
     
    18091810    generate new points and vary their temperature. Temperatures are
    18101811    only updated with iterations in the outer loop. The inner loop is
    1811     over loop over xrange(dwell), and new points are generated for
     1812    over xrange(dwell), and new points are generated for
    18121813    every iteration in the inner loop. (Though whether the proposed
    18131814    new points are accepted is probabilistic.)
     
    18541855    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
    18551856                  learn_rate=learn_rate, lower=lower, upper=upper,
    1856                   m=m, n=n, quench=quench, dwell=dwell)
     1857                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
    18571858
    18581859    current_state, last_state, best_state = _state(), _state(), _state()
     
    18861887                    best_state.x = last_state.x.copy()
    18871888                    best_state.cost = last_state.cost
     1889        if dlg:
     1890            GoOn = dlg.Update(best_state.cost*100,
     1891                newmsg='%s%8.3f\n%s%8.3f%s'%('Temperature =',schedule.T,'MC/SA Residual =',best_state.cost*100,'%'))[0]
     1892            if not GoOn:
     1893                break
    18881894        schedule.update_temp()
    18891895        iters += 1
     
    19391945            lower.append(limits[0])
    19401946            upper.append(limits[1])
    1941            
     1947                       
    19421948    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
    19431949        parmDict[pfx+'Atype'] = item['atType']
     
    19561962        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
    19571963           
    1958     def getRBparms(item,mfx,aTypes,RBdata,atNo,parmDict,varyList):
     1964    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
    19591965        parmDict[mfx+'MolCent'] = item['MolCent']
    19601966        parmDict[mfx+'RBId'] = item['RBId']
     
    19952001        atNo += len(atypes)
    19962002        return atNo
    1997 
    1998     def GetAtomTMX(pfx,RBdata,parmDict):
     2003       
     2004    def GetAtomM(Xdata,SGData):
     2005        Mdata = []
     2006        for xyz in Xdata.T:
     2007            Mdata.append(len(G2spc.GenAtom(xyz,SGData)))
     2008        return np.array(Mdata)
     2009
     2010    def GetAtomTX(RBdata,parmDict):
    19992011        'Needs a doc string'
    2000         atNo = parmDIct['atNo']
     2012        Bmat = parmDict['Bmat']
     2013        atNo = parmDict['atNo']
    20012014        nfixAt = parmDict['nfixAt']
    20022015        Tdata = atNo*[' ',]
    2003         Mdata = np.zeros(atNo)
    2004         Fdata = np.zeros(atNo)
    20052016        Xdata = np.zeros((3,atNo))
    2006         keys = {'Atype:':Tdata,'Amul:':Mdata,
    2007             'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2]}
    2008         nObjs = parmDict['nObjs']
     2017        keys = {':Atype':Tdata,':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
    20092018        for iatm in range(nfixAt):
    20102019            for key in keys:
    2011                 parm = pfx+key+str(iatm)
     2020                parm = ':'+str(iatm)+key
    20122021                if parm in parmDict:
    20132022                    keys[key][iatm] = parmDict[parm]
    2014         return Tdata,Mdata,Xdata
    2015    
    2016 #    def mcsaSfCalc(refList,RBdata,G,SGData,parmDict):
     2023        iatm = nfixAt
     2024        for iObj in range(parmDict['nObj']):
     2025            pfx = str(iObj)+':'
     2026            if parmDict[pfx+'Type'] in ['Vector','Residue']:
     2027                if parmDict[pfx+'Type'] == 'Vector':
     2028                    RBId = parmDict[pfx+':RBId']
     2029                    RBRes = RBdata['Vector'][RBId]
     2030                    aTypes = RBRes['rbTypes']
     2031                    vecs = RBRes['rbVect']
     2032                    mags = RBRes['VectMag']
     2033                    Cart = np.zeros_like(vecs[0])
     2034                    for vec,mag in zip(vecs,mags):
     2035                        Cart += vec*mag
     2036                elif parmDict[pfx+'Type'] == 'Residue':
     2037                    RBId = parmDict[pfx+':RBId']
     2038                    RBRes = RBdata['Residue'][RBId]
     2039                    aTypes = RBRes['rbTypes']
     2040                    Cart = np.array(RBRes['rbXYZ'])
     2041                    for itor,seq in enumerate(RBRes['rbSeq']):
     2042                        tName = pfx+':Tor'+str(itor)
     2043                        QuatA = AVdeg2Q(parmDict[tName],Cart[seq[0]]-Cart[seq[1]])
     2044                        for ride in seq[3]:
     2045                            Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
     2046                if parmDict[pfx+':MolCent'][1]:
     2047                    Cart -= parmDict[pfx+':MolCent'][0]
     2048                Qori = np.array([parmDict[pfx+':Qa'],parmDict[pfx+':Qi'],parmDict[pfx+':Qj'],parmDict[pfx+':Qk']])
     2049                Pos = np.array([parmDict[pfx+':Px'],parmDict[pfx+':Py'],parmDict[pfx+':Pz']])
     2050                for i,x in enumerate(Cart):
     2051                    X = np.inner(Bmat,prodQVQ(Qori,x))+Pos
     2052                    for j in range(3):
     2053                        Xdata[j][iatm] = X[j]
     2054                    Tdata[iatm] = aTypes[i]
     2055                    iatm += 1
     2056            elif parmDict[pfx+'Type'] == 'Atom':
     2057                atNo = parmDict[pfx+'atNo']
     2058                afx = pfx+str(atNo)
     2059                for key in keys:
     2060                    parm = afx+key
     2061                    if parm in parmDict:
     2062                        keys[key][atNo] = parmDict[parm]
     2063            else:
     2064                continue        #skips March Dollase
     2065        return Tdata,Xdata
     2066   
     2067    def calcMDcorr(MDval,MDaxis,Uniq,G):
     2068        sumMD = 0
     2069        for H in Uniq:           
     2070            cosP,sinP = G2lat.CosSinAngle(H,MDaxis,G)
     2071            A = 1.0/np.sqrt((MDval*cosP)**2+sinP**2/MDval)
     2072            sumMD += A**3
     2073        return sumMD
     2074       
     2075    def mcsaCalc(values,refList,rcov,ifInv,RBdata,varyList,parmDict):
    20172076#        ''' Compute structure factors for all h,k,l for phase
    20182077#        input:
    20192078#            refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
    2020 #            G:      reciprocal metric tensor
    2021 #            SGData: space group info. dictionary output from SpcGroup
    20222079#            ParmDict:
    20232080#        puts result F^2 in each ref[8] in refList
    20242081#        '''       
    2025 #        twopi = 2.0*np.pi
    2026 #        twopisq = 2.0*np.pi**2
    2027 #        ast = np.sqrt(np.diag(G))
    2028 #        Mast = twopisq*np.multiply.outer(ast,ast)
    2029 #        Tdata,Mdata,Fdata,Xdata = GetAtomFX(pfx,calcControls,parmDict)
    2030 #        FF = np.zeros(len(Tdata))
    2031 #        if 'N' in parmDict[hfx+'Type']:
    2032 #            FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
    2033 #        else:
    2034 #            FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
    2035 #            FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    2036 #        Uij = np.array(G2lat.U6toUij(Uijdata))
    2037 #        bij = Mast*Uij.T
    2038 #        for refl in refList:
    2039 #            fbs = np.array([0,0])
    2040 #            H = refl[:3]
    2041 #            SQ = 1./(2.*refl[4])**2
    2042 #            SQfactor = 4.0*SQ*twopisq
    2043 #            Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
    2044 #            if not len(refl[-1]):                #no form factors
    2045 #                if 'N' in parmDict[hfx+'Type']:
    2046 #                    refl[-1] = getBLvalues(BLtables)
    2047 #                else:       #'X'
    2048 #                    refl[-1] = getFFvalues(FFtables,SQ)
    2049 #            for i,El in enumerate(Tdata):
    2050 #                FF[i] = refl[-1][El]           
    2051 #            Uniq = refl[11]
    2052 #            phi = refl[12]
    2053 #            phase = twopi*(np.inner(Uniq,(Xdata.T))+phi[:,np.newaxis])
    2054 #            sinp = np.sin(phase)
    2055 #            cosp = np.cos(phase)
    2056 #            occ = Mdata*Fdata/len(Uniq)
    2057 #            fa = np.array([(FF+FP-Bab)*occ*cosp,-FPP*occ*sinp])
    2058 #            fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
    2059 #            if not SGData['SGInv']:
    2060 #                fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
    2061 #                fbs = np.sum(np.sum(fb,axis=1),axis=1)
    2062 #            fasq = fas**2
    2063 #            fbsq = fbs**2        #imaginary
    2064 #            refl[9] = np.sum(fasq)+np.sum(fbsq)
    2065 #            refl[10] = atan2d(fbs[0],fas[0])
     2082        twopi = 2.0*np.pi
     2083        parmDict.update(dict(zip(varyList,values)))
     2084        Tdata,Xdata = GetAtomTX(RBdata,parmDict)
     2085        Mdata = parmDict['Mdata']
     2086        FF = np.zeros(len(Tdata))
     2087        MDval = parmDict['0:MDval']
     2088        MDaxis = parmDict['0:MDaxis']
     2089        Gmat = parmDict['Gmat']
     2090        Srefs = np.zeros(len(refList))
     2091        sumFcsq = 0
     2092        for refl in refList:
     2093            fbs = 0
     2094            H = refl[:3]
     2095            for i,El in enumerate(Tdata):
     2096                FF[i] = refl[7][El]           
     2097            Uniq = refl[8]
     2098            phi = refl[9]
     2099            phase = twopi*(np.inner(Uniq,(Xdata.T))+phi[:,np.newaxis])
     2100            sinp = np.sin(phase)
     2101            cosp = np.cos(phase)
     2102            occ = Mdata/len(Uniq)
     2103            fa = np.asarray(FF*occ*cosp)
     2104            fas = np.sum(fa)
     2105            if not ifInv:
     2106                fb = np.asarray(FF*occ*sinp)
     2107                fbs = np.sum(fb)
     2108            fcsq = (fas**2+fbs**2)*refl[3]*calcMDcorr(MDval,MDaxis,Uniq,Gmat)
     2109            sumFcsq += fcsq
     2110            refl[5] = fcsq
     2111        scale = (parmDict['sumFosq']/sumFcsq)
     2112        for iref,refl in enumerate(refList):
     2113            refl[5] *= scale
     2114            Srefs[iref] = refl[4]-refl[5]
     2115        M = np.inner(Srefs,np.inner(rcov,Srefs))
     2116        return M/parmDict['sumFosq']**2
    20662117   
    20672118    sq8ln2 = np.sqrt(8*np.log(2))
     
    20692120    sq4pi = np.sqrt(4*np.pi)
    20702121    generalData = data['General']
     2122    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     2123    Gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])[0]
    20712124    SGData = generalData['SGData']
    20722125    fixAtoms = data['Atoms']                       #if any
    20732126    cx,ct,cs = generalData['AtomPtrs'][:3]
    20742127    aTypes = set([])
    2075     parmDict = {}
     2128    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
    20762129    varyList = []
    20772130    atNo = 0
     
    20872140        atNo += 1
    20882141    parmDict['nfixAt'] = len(fixAtoms)       
    2089     mcsaControls = generalData['MCSA controls']
    2090     reflName = mcsaControls['Data source']
     2142    MCSA = generalData['MCSA controls']
     2143    reflName = MCSA['Data source']
    20912144    phaseName = generalData['Name']
    20922145    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
     
    21012154            pfx = mfx+str(atNo)+':'
    21022155            getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList)
     2156            parmDict[mfx+'atNo'] = atNo
    21032157            atNo += 1
    21042158        elif item['Type'] in ['Residue','Vector']:
    21052159            pfx = mfx+':'
    2106             atNo = getRBparms(item,pfx,aTypes,RBdata,atNo,parmDict,varyList)
     2160            atNo = getRBparms(item,pfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
    21072161    parmDict['atNo'] = atNo                 #total no. of atoms
    21082162    parmDict['nObj'] = len(MCSAObjs)
     2163    Tdata,Xdata = GetAtomTX(RBdata,parmDict)
     2164    parmDict['Mdata'] = GetAtomM(Xdata,SGData)
    21092165    FFtables = G2el.GetFFtable(aTypes)
    21102166    refs = []
     2167    sumFosq = 0
    21112168    if 'PWDR' in reflName:
    21122169        for ref in reflData:
    21132170            h,k,l,m,d,pos,sig,gam,f = ref[:9]
    2114             if d >= mcsaControls['dmin']:
     2171            if d >= MCSA['dmin']:
    21152172                sig = gamFW(sig,gam)/sq8ln2        #--> sig from FWHM
    21162173                SQ = 0.25/d**2
     
    21182175                FFs = G2el.getFFvalues(FFtables,SQ)
    21192176                refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi])
     2177                sumFosq += f*m
    21202178        nRef = len(refs)
    21212179        rcov = np.zeros((nRef,nRef))
     
    21382196        for iref,refI in enumerate(reflData):
    21392197            h,k,l,m,d,v,f,s = refI
    2140             if d >= mcsaControls['dmin'] and v:       #skip unrefined ones
     2198            if d >= MCSA['dmin'] and v:       #skip unrefined ones
    21412199                SQ = 0.25/d**2
    21422200                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
    21432201                FFs = G2el.getFFvalues(FFtables,SQ)
    21442202                refs.append([h,k,l,m,f*m,iref,0.,FFs,Uniq,phi])
     2203                sumFosq += f*m
    21452204        nRef = len(refs)
    21462205        pfx = str(data['pId'])+'::PWLref:'
     
    21692228        for ref in reflData:
    21702229            [h,k,l,m,d],f = ref[:5],ref[6]
    2171             if d >= mcsaControls['dmin']:
     2230            if d >= MCSA['dmin']:
    21722231                SQ = 0.25/d**2
    21732232                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
    21742233                FFs = G2el.getFFvalues(FFtables,SQ)
    21752234                refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi])
     2235                sumFosq += f*m
    21762236        rcov = np.identity(len(refs))
    2177        
    2178     for parm in parmDict:
    2179         print parm,parmDict[parm]
     2237    parmDict['sumFosq'] = sumFosq
     2238    x0 = [parmDict[val] for val in varyList]
     2239    ifInv = SGData['SGInv']
     2240    results = anneal(mcsaCalc,x0,args=(refs,rcov,ifInv,RBdata,varyList,parmDict),
     2241        schedule=MCSA['Algorithm'], full_output=True,maxiter=MCSA['nRuns'],
     2242        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
     2243        boltzmann=MCSA['boltzmann'], learn_rate=0.5, feps=MCSA['Annealing'][3],
     2244        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
     2245        lower=lower, upper=upper, slope=MCSA['log slope'],dlg=pgbar)
     2246    print results
     2247           
     2248#    parmDict.update(zip(varylist,results[0]))           
    21802249                       
    2181 #    XYZ,aTypes = UpdateMCSAxyz(Bmat,MCSA)       
    2182            
    2183 #    generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50,1.e-6],
    2184 #    'dmin':2.0,'Algorithm':'fast','Jump coeff':[0.95,0.5],'nRuns':1,'boltzmann':1.0,
    2185 #    'fast parms':[1.0,1.0,1.0],'log slope':0.9}
    2186 
    21872250    return {}
    21882251
  • trunk/GSASIIphsGUI.py

    r946 r949  
    127127            'MCSA controls' not in generalData:
    128128            generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50,1.e-6],
    129             'dmin':2.0,'Algorithm':'fast','Jump coeff':[0.95,0.5],'nRuns':1,'boltzmann':1.0,
     129            'dmin':2.0,'Algorithm':'fast','Jump coeff':[0.95,0.5],'nRuns':50,'boltzmann':1.0,
    130130            'fast parms':[1.0,1.0,1.0],'log slope':0.9}
    131131# end of patches
     
    726726            mcsaSizer.Add((5,5),)
    727727            line2Sizer = wx.BoxSizer(wx.HORIZONTAL)
    728             Rchoice = ['1','2','3','5','10','15','20','50','100','200','500']
    729             line2Sizer.Add(wx.StaticText(General,label=' No. MC/SA runs: '),0,wx.ALIGN_CENTER_VERTICAL)
     728            Rchoice = ['10','15','20','50','100','200','500']
     729            line2Sizer.Add(wx.StaticText(General,label=' No. temp. steps: '),0,wx.ALIGN_CENTER_VERTICAL)
    730730            noRuns = wx.ComboBox(General,-1,value=str(MCSA.get('nRuns',1)),choices=Rchoice,
    731731                style=wx.CB_READONLY|wx.CB_DROPDOWN)
Note: See TracChangeset for help on using the changeset viewer.