Changeset 951


Ignore:
Timestamp:
Jun 14, 2013 2:38:20 PM (8 years ago)
Author:
vondreele
Message:

more on MC/SA

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIgrid.py

    r939 r951  
    8080] = [wx.NewId() for item in range(4)]
    8181
    82 [ wxID_ADDMCSAATOM,wxID_ADDMCSARB,wxID_CLEARMCSARB,wxID_MOVEMCSA,
    83 ] = [wx.NewId() for item in range(4)]
     82[ wxID_ADDMCSAATOM,wxID_ADDMCSARB,wxID_CLEARMCSARB,wxID_MOVEMCSA,wxID_MCSACLEARRESULTS,
     83] = [wx.NewId() for item in range(5)]
    8484
    8585[ wxID_CLEARTEXTURE,wxID_REFINETEXTURE,
     
    17421742            text='Clear map')
    17431743        self.GeneralCalc.Append(help='Run Monte Carlo - Simulated Annealing',id=wxID_RUNMCSA, kind=wx.ITEM_NORMAL,
    1744             text='Run MC/SA')
     1744            text='MC/SA')
    17451745        self.PostfillDataMenu()
    17461746       
     
    18731873        self.MCSAEdit.Append(id=wxID_MOVEMCSA, kind=wx.ITEM_NORMAL,text='Move MC/SA solution',
    18741874            help='Move MC/SA solution to atom list' )
     1875        self.MCSAEdit.Append(id=wxID_MCSACLEARRESULTS, kind=wx.ITEM_NORMAL,text='Clear results',
     1876            help='Clear table of MC/SA results' )
    18751877        self.PostfillDataMenu()
    18761878           
  • trunk/GSASIImath.py

    r950 r951  
    24602460            if parmDict[pfx+'Type'] in ['Vector','Residue']:
    24612461                if parmDict[pfx+'Type'] == 'Vector':
    2462                     RBId = parmDict[pfx+':RBId']
     2462                    RBId = parmDict[pfx+'RBId']
    24632463                    RBRes = RBdata['Vector'][RBId]
    24642464                    aTypes = RBRes['rbTypes']
     
    24692469                        Cart += vec*mag
    24702470                elif parmDict[pfx+'Type'] == 'Residue':
    2471                     RBId = parmDict[pfx+':RBId']
     2471                    RBId = parmDict[pfx+'RBId']
    24722472                    RBRes = RBdata['Residue'][RBId]
    24732473                    aTypes = RBRes['rbTypes']
    24742474                    Cart = np.array(RBRes['rbXYZ'])
    24752475                    for itor,seq in enumerate(RBRes['rbSeq']):
    2476                         tName = pfx+':Tor'+str(itor)
     2476                        tName = pfx+'Tor'+str(itor)
    24772477                        QuatA = AVdeg2Q(parmDict[tName],Cart[seq[0]]-Cart[seq[1]])
    24782478                        for ride in seq[3]:
    24792479                            Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
    2480                 if parmDict[pfx+':MolCent'][1]:
    2481                     Cart -= parmDict[pfx+':MolCent'][0]
    2482                 Qori = np.array([parmDict[pfx+':Qa'],parmDict[pfx+':Qi'],parmDict[pfx+':Qj'],parmDict[pfx+':Qk']])
    2483                 Pos = np.array([parmDict[pfx+':Px'],parmDict[pfx+':Py'],parmDict[pfx+':Pz']])
     2480                if parmDict[pfx+'MolCent'][1]:
     2481                    Cart -= parmDict[pfx+'MolCent'][0]
     2482                Qori = np.array([parmDict[pfx+'Qa'],parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
     2483                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
    24842484                for i,x in enumerate(Cart):
    24852485                    X = np.inner(Bmat,prodQVQ(Qori,x))+Pos
     
    24902490            elif parmDict[pfx+'Type'] == 'Atom':
    24912491                atNo = parmDict[pfx+'atNo']
    2492                 afx = pfx+str(atNo)
    24932492                for key in keys:
    2494                     parm = afx+key
     2493                    parm = pfx+key[1:]              #remove extra ':'
    24952494                    if parm in parmDict:
    24962495                        keys[key][atNo] = parmDict[parm]
     
    25752574    parmDict['nfixAt'] = len(fixAtoms)       
    25762575    MCSA = generalData['MCSA controls']
    2577     Results = MCSA.get('Results',[])
    25782576    reflName = MCSA['Data source']
    25792577    phaseName = generalData['Name']
     
    25872585            getMDparms(item,mfx,parmDict,varyList)
    25882586        elif item['Type'] == 'Atom':
    2589             pfx = mfx+str(atNo)+':'
    2590             getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList)
     2587            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
    25912588            parmDict[mfx+'atNo'] = atNo
    25922589            atNo += 1
    25932590        elif item['Type'] in ['Residue','Vector']:
    2594             pfx = mfx+':'
    2595             atNo = getRBparms(item,pfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
     2591            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
    25962592    parmDict['atNo'] = atNo                 #total no. of atoms
    25972593    parmDict['nObj'] = len(MCSAObjs)
     
    26732669    x0 = [parmDict[val] for val in varyList]
    26742670    ifInv = SGData['SGInv']
    2675     for i in range(MCSA['Cycles']):     
    2676         results = anneal(mcsaCalc,x0,args=(refs,rcov,ifInv,RBdata,varyList,parmDict),
    2677             schedule=MCSA['Algorithm'], full_output=True,maxiter=MCSA['nRuns'],
    2678             T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
    2679             boltzmann=MCSA['boltzmann'], learn_rate=0.5, feps=MCSA['Annealing'][3],
    2680             quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
    2681             lower=lower, upper=upper, slope=MCSA['log slope'],dlg=pgbar)
    2682         Results.append([results[1],results[2],results[0],varyList])
     2671    results = anneal(mcsaCalc,x0,args=(refs,rcov,ifInv,RBdata,varyList,parmDict),
     2672        schedule=MCSA['Algorithm'], full_output=True,
     2673        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
     2674        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
     2675        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
     2676        lower=lower, upper=upper, slope=MCSA['log slope'],dlg=pgbar)
     2677    return [False,results[1],results[2],results[0],varyList]
    26832678
    26842679       
  • trunk/GSASIIphsGUI.py

    r950 r951  
    7070
    7171    :param wx.frame G2frame: the main GSAS-II frame object
    72 
    7372    :param wx.TreeItemId Item: the tree item that was selected
    74 
    7573    :param dict data: all the information on the phase in a dictionary
    76 
    7774    :param int oldPage: This sets a tab to select when moving
    7875      from one phase to another, in which case the same tab is selected
     
    8683    if 'RBModels' not in data:
    8784        data['RBModels'] = {}
     85    if isinstance(data['MCSA']['Results'],dict):
     86        data['MCSA']['Results'] = []
    8887    if 'MCSA' not in data:
    89         data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[0.,3.],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}}
     88        data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[0.3,3.],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}}
    9089#end patch   
    9190
     
    126125        if 'Algolrithm' in generalData.get('MCSA controls',{}) or \
    127126            'MCSA controls' not in generalData:
    128             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':50,'boltzmann':1.0,
     127            generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50],
     128            'dmin':2.0,'Algorithm':'fast','Jump coeff':[0.95,0.5],'boltzmann':1.0,
    130129            'fast parms':[1.0,1.0,1.0],'log slope':0.9,'Cycles':1,'Results':[]}
    131130# end of patches
     
    659658
    660659            def OnCycles(event):
    661                 MCSA['Cycles'] = int(noRuns.GetValue())
     660                MCSA['Cycles'] = int(cycles.GetValue())
    662661                               
    663             def OnNoRuns(event):
    664                 MCSA['nRuns'] = int(noRuns.GetValue())
    665            
    666662            def OnAlist(event):
    667663                MCSA['Algorithm'] = Alist.GetValue()
     
    697693                            MCSA['Annealing'][ind] = val
    698694                    except ValueError:
    699                         pass
     695                        Obj.SetValue(fmt%(MCSA['Annealing'][ind]))
    700696                else:
    701697                    try:
     
    703699                        if .0 <= val:
    704700                            MCSA['Annealing'][ind] = val
     701                        Obj.SetValue(fmt%(MCSA['Annealing'][ind]))
    705702                    except ValueError:
    706                         pass                   
    707                 Obj.SetValue(fmt%(MCSA['Annealing'][ind]))
     703                        MCSA['Annealing'][ind] = None                   
     704                        Obj.SetValue(str(MCSA['Annealing'][ind]))
    708705                       
    709 #            MCSA = {'boltzmann':1.0,
    710706            refList = []
    711707            if len(data['Pawley ref']):
     
    726722            dmin.Bind(wx.EVT_KILL_FOCUS,OnDmin)
    727723            lineSizer.Add(dmin,0,wx.ALIGN_CENTER_VERTICAL)
    728             lineSizer.Add(wx.StaticText(General,label=' Cycles: '),0,wx.ALIGN_CENTER_VERTICAL)
     724            mcsaSizer.Add(lineSizer)
     725            mcsaSizer.Add((5,5),)
     726            line2Sizer = wx.BoxSizer(wx.HORIZONTAL)
     727            line2Sizer.Add(wx.StaticText(General,label=' Cycles: '),0,wx.ALIGN_CENTER_VERTICAL)
    729728            Cchoice = ['1','2','3','5','10','15','20','30']
    730729            cycles = wx.ComboBox(General,-1,value=str(MCSA.get('Cycles',1)),choices=Cchoice,
    731730                style=wx.CB_READONLY|wx.CB_DROPDOWN)
    732731            cycles.Bind(wx.EVT_COMBOBOX,OnCycles)       
    733             lineSizer.Add(cycles,0,wx.ALIGN_CENTER_VERTICAL)
    734             mcsaSizer.Add(lineSizer)
    735             mcsaSizer.Add((5,5),)
    736             line2Sizer = wx.BoxSizer(wx.HORIZONTAL)
    737             Rchoice = ['10','15','20','50','100','200','500']
    738             line2Sizer.Add(wx.StaticText(General,label=' No. temp. steps: '),0,wx.ALIGN_CENTER_VERTICAL)
    739             noRuns = wx.ComboBox(General,-1,value=str(MCSA.get('nRuns',1)),choices=Rchoice,
    740                 style=wx.CB_READONLY|wx.CB_DROPDOWN)
    741             noRuns.Bind(wx.EVT_COMBOBOX,OnNoRuns)
    742             line2Sizer.Add(noRuns,0,wx.ALIGN_CENTER_VERTICAL)
     732            line2Sizer.Add(cycles,0,wx.ALIGN_CENTER_VERTICAL)
    743733            Achoice = ['log','fast','cauchy','boltzmann','Tremayne']
    744734            line2Sizer.Add(wx.StaticText(General,label=' MC/SA schedule: '),0,wx.ALIGN_CENTER_VERTICAL)
     
    772762            line3Sizer = wx.BoxSizer(wx.HORIZONTAL)
    773763            line3Sizer.Add(wx.StaticText(General,label=' Annealing schedule: '),0,wx.ALIGN_CENTER_VERTICAL)
    774             names = [' Start temp: ',' Final temp: ',' No. trials: ',' feps: ']
    775             fmts = ['%.1f','%.5f','%d','%.2g']
     764            names = [' Start temp: ',' Final temp: ',' No. trials: ']
     765            fmts = ['%.1f','%.5f','%d']
    776766            for i,[name,fmt] in enumerate(zip(names,fmts)):
     767                if MCSA['Annealing'][i]:
     768                    text = fmt%(MCSA['Annealing'][i])
     769                else:
     770                    text = 'None'
    777771                line3Sizer.Add(wx.StaticText(General,label=name),0,wx.ALIGN_CENTER_VERTICAL)
    778                 anneal =  wx.TextCtrl(General,-1,value=fmt%(MCSA['Annealing'][i]),style=wx.TE_PROCESS_ENTER)
     772                anneal =  wx.TextCtrl(General,-1,value=text,style=wx.TE_PROCESS_ENTER)
    779773                anneal.Bind(wx.EVT_TEXT_ENTER,OnAnneal)       
    780774                anneal.Bind(wx.EVT_KILL_FOCUS,OnAnneal)
     
    41614155                try:
    41624156                    rmin,rmax = [float(Range[i]) for i in range(2)]
    4163                     if rmin >= rmax:
     4157                    if 0. < rmin < rmax:
     4158                        pass
     4159                    else:
    41644160                        raise ValueError
    41654161                except (ValueError,IndexError):
     
    42024198            poSizer.Add(poAxis,0,wx.ALIGN_CENTER_VERTICAL)
    42034199            return poSizer
     4200           
     4201        def ResultsSizer(Results):
     4202           
     4203            def OnCellChange(event):
     4204                r,c = event.GetRow(),event.GetCol()
     4205                if c == 0:
     4206                    Models = data['MCSA']['Models']
     4207                    for row in range(resultsGrid.GetNumberRows()):
     4208                        resultsTable.SetValue(row,c,False)
     4209                    resultsTable.SetValue(r,c,True)
     4210                    resultsGrid.ForceRefresh()
     4211                    result = Results[r]
     4212                    for key,val in zip(result[4],result[3]):
     4213                        vals = key.split(':')
     4214                        nObj,name = int(vals[0]),vals[1]
     4215                        if 'A' in name:
     4216                            ind = ['Ax','Ay','Az'].index(name)
     4217                            Models[nObj]['Pos'][0][ind] = val                           
     4218                        elif 'O' in name:
     4219                            ind = ['Oa','Oi','Oj','Ok'].index(name)
     4220                            Models[nObj]['Ori'][0][ind] = val                           
     4221                        elif 'P' in name:
     4222                            ind = ['Px','Py','Pz'].index(name)
     4223                            Models[nObj]['Pos'][0][ind] = val                           
     4224                        elif 'T' in name:
     4225                            tnum = int(name.split('Tor'))
     4226                            Models[nObj]['Tor'][0][tnum] = val                                                       
     4227                        else:       #March Dollase
     4228                            Models[0]['Coef'][0] = val
     4229                    UpdateMCSA()
     4230                    G2plt.PlotStructure(G2frame,data)
     4231               
     4232            resultsSizer = wx.BoxSizer(wx.VERTICAL)
     4233            maxVary = 0
     4234            resultVals = []
     4235            for result in Results:
     4236                maxVary = max(maxVary,len(result[3]))
     4237                resultVals.append(result[:3]+list(result[3]))
     4238            rowLabels = []
     4239            for i in range(len(Results)): rowLabels.append(str(i))
     4240            colLabels = ['Select','Residual','Tmin',]
     4241            for i in range(maxVary): colLabels.append('variable:'+str(i))
     4242            Types = [wg.GRID_VALUE_BOOL,wg.GRID_VALUE_FLOAT+':10,4',
     4243                wg.GRID_VALUE_FLOAT+':10,4',]+maxVary*[wg.GRID_VALUE_FLOAT+':10,5',]
     4244            resultsTable = G2gd.Table(resultVals,rowLabels=rowLabels,colLabels=colLabels,types=Types)
     4245            resultsGrid = G2gd.GSGrid(MCSA)
     4246            resultsGrid.SetTable(resultsTable, True)
     4247            resultsGrid.Bind(wg.EVT_GRID_CELL_LEFT_CLICK, OnCellChange)
     4248            resultsGrid.AutoSizeColumns(True)
     4249            for r in range(resultsGrid.GetNumberRows()):
     4250                for c in range(resultsGrid.GetNumberCols()):
     4251                    if c == 0:
     4252                        resultsGrid.SetReadOnly(r,c,isReadOnly=False)
     4253                    else:
     4254                        resultsGrid.SetCellStyle(r,c,VERY_LIGHT_GREY,True)
     4255            resultsSizer.Add(resultsGrid)
     4256            return resultsSizer
     4257           
    42044258       
    42054259        # UpdateMCSA executable code starts here
     
    42384292            mainSizer.Add((5,5),0)
    42394293        else:
    4240             for result in data['MCSA']['Results']:
    4241                 print result
     4294            mainSizer.Add((5,5),0)
     4295            mainSizer.Add(wx.StaticText(MCSA,-1,'MC/SA results:'),0,wx.ALIGN_CENTER_VERTICAL)
     4296            mainSizer.Add((5,5),0)
     4297            Results = data['MCSA']['Results']
     4298            mainSizer.Add(ResultsSizer(Results))
    42424299
    42434300        MCSA.SetSizer(mainSizer)
     
    42944351        pgbar.SetSize(Size)
    42954352        try:
    4296             G2mth.mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar)
     4353            for i in range(mcsaControls['Cycles']):
     4354                MCSAdata['Results'].append(G2mth.mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar))
    42974355        finally:
    42984356            pgbar.Destroy()
    4299         if not data['Drawing']:                 #if new drawing - no drawing data!
    4300             SetupDrawingData()
    43014357        UpdateMCSA()
     4358        G2plt.PlotStructure(G2frame,data)
    43024359
    43034360    def OnMCSAaddAtom(event):
     
    43634420            AtomAdd(x,y,z,atype,Name=atype+'(%d)'%(iat+1))           
    43644421        G2plt.PlotStructure(G2frame,data)
     4422       
     4423    def OnClearResults(event):
     4424        data['MCSA']['Results'] = []
     4425        UpdateMCSA()
    43654426                   
    43664427################################################################################
     
    49465007            G2frame.dataFrame.Bind(wx.EVT_MENU, OnMCSAclear, id=G2gd.wxID_CLEARMCSARB)
    49475008            G2frame.dataFrame.Bind(wx.EVT_MENU, OnMCSAmove, id=G2gd.wxID_MOVEMCSA)
     5009            G2frame.dataFrame.Bind(wx.EVT_MENU, OnClearResults, id=G2gd.wxID_MCSACLEARRESULTS)
    49485010            UpdateMCSA()                       
    49495011            wx.CallAfter(G2plt.PlotStructure,G2frame,data)
  • trunk/GSASIIplot.py

    r939 r951  
    26162616    '''
    26172617
     2618    def FindPeaksBonds(XYZ):
     2619        rFact = drawingData['radiusFactor']
     2620        Bonds = [[] for x in XYZ]
     2621        for i,xyz in enumerate(XYZ):
     2622            Dx = XYZ-xyz
     2623            dist = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=1))
     2624            IndB = ma.nonzero(ma.masked_greater(dist,rFact*2.2))
     2625            for j in IndB[0]:
     2626                Bonds[i].append(Dx[j]/2.)
     2627                Bonds[j].append(-Dx[j]/2.)
     2628        return Bonds
     2629
    26182630    ForthirdPI = 4.0*math.pi/3.0
    26192631    generalData = data['General']
     
    26242636    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
    26252637    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
     2638    SGData = generalData['SGData']
    26262639    Mydir = generalData['Mydir']
    26272640    atomData = data['Atoms']
     
    26432656    MCSA = data.get('MCSA',{})
    26442657    mcsaModels = MCSA.get('Models',[])
     2658    if mcsaModels:
     2659            mcsaXYZ,atTypes = G2mth.UpdateMCSAxyz(Bmat,MCSA)
     2660            XYZeq = []
     2661            for xyz in mcsaXYZ:
     2662                XYZeq += G2spc.GenAtom(xyz,SGData)[0][1:]       #skip self xyz
     2663           
     2664            mcsaBonds = FindPeaksBonds(mcsaXYZ)       
    26452665    drawAtoms = drawingData.get('Atoms',[])
    26462666    mapData = {}
     
    26732693    ctrlDown = False
    26742694   
    2675     def FindPeaksBonds(XYZ):
    2676         rFact = drawingData['radiusFactor']
    2677         Bonds = [[] for x in XYZ]
    2678         for i,xyz in enumerate(XYZ):
    2679             Dx = XYZ-xyz
    2680             dist = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=1))
    2681             IndB = ma.nonzero(ma.masked_greater(dist,rFact*2.2))
    2682             for j in IndB[0]:
    2683                 Bonds[i].append(Dx[j]/2.)
    2684                 Bonds[j].append(-Dx[j]/2.)
    2685         return Bonds
    2686 
    26872695    def OnKeyBox(event):
    26882696        import Image
     
    34963504                RenderLabel(x,y,z,name,0.2,Or)
    34973505        if len(mcsaModels) > 1 and pageName == 'MC/SA':             #skip the default MD entry
    3498             XYZ,atTypes = G2mth.UpdateMCSAxyz(Bmat,MCSA)
    3499             rbBonds = FindPeaksBonds(XYZ)
    3500             for ind,[x,y,z] in enumerate(XYZ):
     3506            for ind,[x,y,z] in enumerate(mcsaXYZ):
    35013507                aType = atTypes[ind]
    35023508                name = '  '+aType+str(ind)
    35033509                color = np.array(MCSA['AtInfo'][aType][1])
    35043510                RenderSphere(x,y,z,0.2,color/255.)
    3505                 RenderBonds(x,y,z,rbBonds[ind],0.03,Gr)
     3511                RenderBonds(x,y,z,mcsaBonds[ind],0.03,Gr)
    35063512                RenderLabel(x,y,z,name,0.2,Or)
    35073513        if Backbones:
Note: See TracChangeset for help on using the changeset viewer.