Changeset 399


Ignore:
Timestamp:
Oct 25, 2011 4:53:33 PM (10 years ago)
Author:
vondreele
Message:

Implement sequential refinement
remove print "load" & "save" for each item in Tree
revise application of azimuth offset - azimuths are now all "true" with correction

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASII.py

    r397 r399  
    6464
    6565[wxID_FILECLOSE, wxID_FILEEXIT, wxID_FILEOPEN,  wxID_FILESAVE, wxID_FILESAVEAS,
    66 wxID_REFINE, wxID_SOLVE, wxID_MAKEPDFS, wxID_VIEWLSPARMS,
    67 ] = [wx.NewId() for _init_coll_File_Items in range(9)]
     66wxID_REFINE, wxID_SOLVE, wxID_MAKEPDFS, wxID_VIEWLSPARMS, wxID_SEQREFINE,
     67] = [wx.NewId() for _init_coll_File_Items in range(10)]
    6868
    6969[wxID_PWDRREAD,wxID_SNGLREAD,wxID_ADDPHASE,wxID_DELETEPHASE,
     
    153153        self.Refine.Enable(False)
    154154        self.Bind(wx.EVT_MENU, self.OnRefine, id=wxID_REFINE)
     155        self.SeqRefine = parent.Append(help='', id=wxID_SEQREFINE, kind=wx.ITEM_NORMAL,
     156            text='Sequental refine')
     157        self.SeqRefine.Enable(False)
     158        self.Bind(wx.EVT_MENU, self.OnSeqRefine, id=wxID_SEQREFINE)
    155159        self.Solve = parent.Append(help='', id=wxID_SOLVE, kind=wx.ITEM_NORMAL,
    156160            text='Solve')
     
    299303            self.PatternTree.Expand(self.root)
    300304            self.Refine.Enable(True)
     305            self.SeqRefine.Enable(True)
    301306            self.Solve.Enable(True)
    302307
     
    10771082                            if data:
    10781083                                self.Refine.Enable(True)
     1084                                self.SeqRefine.Enable(True)
    10791085                                self.Solve.Enable(True)         #not right but something needed here
    10801086                        item, cookie = self.PatternTree.GetNextChild(self.root, cookie)               
     
    15171523        finally:
    15181524            dlg.Destroy()
     1525
     1526    def OnSeqRefine(self,event):
     1527        self.OnFileSave(event)
     1528        Id = G2gd.GetPatternTreeItemId(self,self.root,'Sequental results')
     1529        if not Id:
     1530            Id = self.PatternTree.AppendItem(self.root,text='Sequental results')
     1531            self.PatternTree.SetItemPyData(Id,{})           
     1532        dlg = wx.ProgressDialog('Residual for histogram 0','Powder profile Rwp =',101.0,
     1533            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_CAN_ABORT)
     1534        screenSize = wx.ClientDisplayRect()
     1535        Size = dlg.GetSize()
     1536        dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
     1537        try:
     1538            G2str.SeqRefine(self.GSASprojectfile,dlg)
     1539        finally:
     1540            dlg.Destroy()       
     1541        dlg = wx.MessageDialog(self,'Load new result?','Refinement results',wx.OK|wx.CANCEL)
     1542        try:
     1543            if dlg.ShowModal() == wx.ID_OK:
     1544                Id = 0
     1545                self.PatternTree.DeleteChildren(self.root)
     1546                if self.HKL: self.HKL = []
     1547                if self.G2plotNB.plotList:
     1548                    self.G2plotNB.clear()
     1549                G2IO.ProjFileOpen(self)
     1550                item, cookie = self.PatternTree.GetFirstChild(self.root)
     1551                while item and not Id:
     1552                    name = self.PatternTree.GetItemText(item)
     1553                    if name[:4] in ['PWDR','HKLF']:
     1554                        Id = item
     1555                    item, cookie = self.PatternTree.GetNextChild(self.root, cookie)               
     1556                if Id:
     1557                    self.PatternTree.SelectItem(Id)
     1558        finally:
     1559            dlg.Destroy()
    15191560       
    15201561    def OnSolve(self,event):
  • trunk/GSASIIIO.py

    r397 r399  
    801801                break
    802802            datum = data[0]
    803             print 'load: ',datum[0]
    804803           
    805804            Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
     
    809808                self.PatternTree.SetItemPyData(Id,datum[1])
    810809            for datus in data[1:]:
    811                 print '    load: ',datus[0]
    812810                sub = self.PatternTree.AppendItem(Id,datus[0])
    813811                self.PatternTree.SetItemPyData(sub,datus[1])
     
    833831                data = []
    834832                name = self.PatternTree.GetItemText(item)
    835                 print 'save: ',name
    836833                data.append([name,self.PatternTree.GetItemPyData(item)])
    837834                item2, cookie2 = self.PatternTree.GetFirstChild(item)
    838835                while item2:
    839836                    name = self.PatternTree.GetItemText(item2)
    840                     print '    save: ',name
    841837                    data.append([name,self.PatternTree.GetItemPyData(item2)])
    842838                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
  • trunk/GSASIIgrid.py

    r397 r399  
    1111import time
    1212import cPickle
     13import numpy as np
    1314import GSASIIpath
    1415import GSASIIplot as G2plt
     
    461462            return self.rowLabels[row]
    462463       
     464    def GetColValues(self, col):
     465        data = []
     466        for row in range(self.GetNumberRows()):
     467            data.append(self.GetValue(row, col))
     468        return data
     469       
    463470    def GetRowValues(self, row):
    464471        data = []
     
    580587        data['shift factor'] = value
    581588        Factr.SetValue('%.5f'%(value))
     589       
     590    def OnSelectData(event):
     591        choices = ['All',]+GetPatternTreeDataNames(self,['PWDR',])
     592        if 'Seq Data' in data:
     593            sel = []
     594            for item in data['Seq Data']:
     595                sel.append(choices.index(item))
     596        names = []
     597        dlg = wx.MultiChoiceDialog(self,'Select data:','Sequential refinement',choices)
     598        dlg.SetSelections(sel)
     599        if dlg.ShowModal() == wx.ID_OK:
     600            sel = dlg.GetSelections()
     601            for i in sel: names.append(choices[i])
     602            if 'All' in names:
     603                names = choices[1:]
     604               
     605        dlg.Destroy()
     606        data['Seq Data'] = names
     607        reverseSel.Enable(True)
     608       
     609    def OnReverse(event):
     610        data['Reverse Seq'] = reverseSel.GetValue()
    582611       
    583612    if self.dataDisplay:
     
    610639    Factr.Bind(wx.EVT_KILL_FOCUS,OnFactor)
    611640    LSSizer.Add(Factr,0,wx.ALIGN_CENTER_VERTICAL)
    612        
    613641    mainSizer.Add(LSSizer)
    614642    mainSizer.Add((5,5),0)
     643   
     644    SeqSizer = wx.BoxSizer(wx.HORIZONTAL)
     645    SeqSizer.Add(wx.StaticText(self.dataDisplay,label=' Sequential Refinement Powder Data: '),0,wx.ALIGN_CENTER_VERTICAL)
     646    selSeqData = wx.Button(self.dataDisplay,-1,label=' Select data')
     647    selSeqData.Bind(wx.EVT_BUTTON,OnSelectData)
     648    SeqSizer.Add(selSeqData,0,wx.ALIGN_CENTER_VERTICAL)
     649    SeqSizer.Add((5,0),0)
     650    reverseSel = wx.CheckBox(self.dataDisplay,-1,label=' Reverse order?')
     651    reverseSel.Bind(wx.EVT_CHECKBOX,OnReverse)
     652    if 'Seq Data' not in data:
     653        reverseSel.Enable(False)
     654    if 'Reverse Seq' in data:
     655        reverseSel.SetValue(data['Reverse Seq'])
     656    SeqSizer.Add(reverseSel,0,wx.ALIGN_CENTER_VERTICAL)
     657   
     658   
     659    mainSizer.Add(SeqSizer)
     660    mainSizer.Add((5,5),0)
     661       
    615662    mainSizer.Add(wx.StaticText(self.dataDisplay,label=' Density Map Controls:'),0,wx.ALIGN_CENTER_VERTICAL)
    616663
     
    633680            self.dataDisplay.AppendText(line+'\n')
    634681           
     682def UpdateSeqResults(self,data):
     683    if not data:
     684        print 'No sequential refinement results'
     685        return
     686    histNames = data['histNames']
     687   
     688    def ColSelect(event):
     689        cols = self.dataDisplay.GetSelectedCols()
     690        plotData = []
     691        plotNames = []
     692        for col in cols:
     693            plotData.append(self.SeqTable.GetColValues(col))
     694            plotNames.append(self.SeqTable.GetColLabelValue(col))
     695        plotData = np.array(plotData)
     696        G2plt.PlotSeq(self,plotData,plotNames)
     697       
     698    if self.dataDisplay:
     699        self.dataDisplay.Destroy()
     700    self.dataFrame.SetMenuBar(self.dataFrame.BlankMenu)
     701    self.dataFrame.SetLabel('Sequental refinement results')
     702    self.dataFrame.CreateStatusBar()
     703    colLabels = data['varyList']
     704    Types = len(data['varyList'])*[wg.GRID_VALUE_FLOAT,]
     705    seqList = [list(data[name]['variables']) for name in histNames]
     706    self.SeqTable = Table(seqList,colLabels=colLabels,rowLabels=histNames,types=Types)
     707    self.dataDisplay = GSGrid(parent=self.dataFrame)
     708    self.dataDisplay.SetTable(self.SeqTable, True)
     709    self.dataDisplay.EnableEditing(False)
     710    self.dataDisplay.Bind(wg.EVT_GRID_LABEL_LEFT_DCLICK, ColSelect)
     711    self.dataDisplay.SetRowLabelSize(8*len(histNames[0]))       #pretty arbitrary 8
     712    self.dataDisplay.SetMargins(0,0)
     713    self.dataDisplay.AutoSizeColumns(True)
     714    self.dataFrame.setSizePosLeft([700,350])
     715               
    635716def UpdateConstraints(self,data):
    636717    Histograms,Phases = self.GetUsedHistogramsAndPhasesfromTree()
     
    9331014    self.dataDisplay.SetSize(mainSizer.Fit(self.dataFrame))
    9341015    self.dataFrame.setSizePosLeft(mainSizer.Fit(self.dataFrame))
    935        
     1016
     1017def GetPatternTreeDataNames(self,dataTypes):
     1018    names = []
     1019    item, cookie = self.PatternTree.GetFirstChild(self.root)       
     1020    while item:
     1021        name = self.PatternTree.GetItemText(item)
     1022        if name[:4] in dataTypes:
     1023            names.append(name)
     1024        item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
     1025    return names
    9361026                         
    9371027def GetPatternTreeItemId(self, parentId, itemText):
     
    9971087                self.PatternTree.SetItemPyData(item,data)                             
    9981088            self.Refine.Enable(True)
     1089            self.SeqRefine.Enable(True)
    9991090            UpdateControls(self,data)
     1091        elif self.PatternTree.GetItemText(item) == 'Sequental results':
     1092            data = self.PatternTree.GetItemPyData(item)
     1093            UpdateSeqResults(self,data)           
    10001094        elif self.PatternTree.GetItemText(item) == 'Covariance':
    10011095            data = self.PatternTree.GetItemPyData(item)
  • trunk/GSASIIimgGUI.py

    r398 r399  
    214214        data['rings'] = []
    215215        data['ellipses'] = []
    216         self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMCLEARCALIB,enable=False)   
    217         self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMRECALIBRATE,enable=False)   
     216#        self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMRECALIBRATE,enable=False)   
    218217        G2plt.PlotExposedImage(self,event=event)
    219218           
    220219    def OnCalibrate(event):       
    221         self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMCLEARCALIB,enable=True)   
    222220        self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMRECALIBRATE,enable=True)   
    223221        self.dataFrame.GetStatusBar().SetStatusText('Select > 4 points on 1st used ring; LB to pick, RB on point to delete else RB to finish')
     
    399397    self.dataFrame.Bind(wx.EVT_MENU, OnClearCalib, id=G2gd.wxID_IMCLEARCALIB)
    400398    if not data['rings']:
    401         self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMCLEARCALIB,enable=False)   
    402399        self.dataFrame.ImageEdit.Enable(id=G2gd.wxID_IMRECALIBRATE,enable=False)   
    403400    self.dataFrame.Bind(wx.EVT_MENU, OnIntegrate, id=G2gd.wxID_IMINTEGRATE)
  • trunk/GSASIIphsGUI.py

    r397 r399  
    20932093            #how about HKLF data? This is only for PWDR data
    20942094            Obj = event.GetEventObject()
    2095             Obj.SetValue(False)
    20962095            hist = Indx[Obj.GetId()]
    20972096            sourceDict = UseList[hist]
     
    23382337            showData.Bind(wx.EVT_CHECKBOX, OnShowData)
    23392338            showSizer.Add(showData,0,wx.ALIGN_CENTER_VERTICAL)
    2340             copyData = wx.CheckBox(dataDisplay,-1,label=' Copy?')
     2339            copyData = wx.Button(dataDisplay,-1,label=' Copy?')
    23412340            Indx[copyData.GetId()] = item
    2342             copyData.Bind(wx.EVT_CHECKBOX,OnCopyData)
     2341            copyData.Bind(wx.EVT_BUTTON,OnCopyData)
    23432342            showSizer.Add(copyData,wx.ALIGN_CENTER_VERTICAL)
    23442343            mainSizer.Add(showSizer,0,wx.ALIGN_CENTER_VERTICAL)
  • trunk/GSASIIplot.py

    r396 r399  
    14331433    Plot.set_ylabel('Variable name')
    14341434    Page.canvas.draw()
    1435 
     1435   
     1436def PlotSeq(self,SeqData,SeqNames):
     1437   
     1438    try:
     1439        plotNum = self.G2plotNB.plotList.index('Sequential refinement')
     1440        Page = self.G2plotNB.nb.GetPage(plotNum)
     1441        Page.figure.clf()
     1442        Plot = Page.figure.gca()
     1443        if not Page.IsShown():
     1444            Page.Show()
     1445    except ValueError:
     1446        Plot = self.G2plotNB.addMpl('Sequential refinement').gca()
     1447        plotNum = self.G2plotNB.plotList.index('Sequential refinement')
     1448        Page = self.G2plotNB.nb.GetPage(plotNum)
     1449       
     1450    Page.SetFocus()
     1451    self.G2plotNB.status.SetFields(['',''])
     1452    if len(SeqData):   
     1453        X = np.arange(0,len(SeqData[0]),1)
     1454        for Y,name in zip(SeqData,SeqNames):
     1455            Plot.plot(X,Y,label=name)       
     1456        Plot.legend(loc='best')
     1457        Page.canvas.draw()
    14361458           
    14371459def PlotExposedImage(self,newPlot=False,event=None):
     
    18111833                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
    18121834            for i in range(Nazm):
    1813                 cake = LRAzim[0]+i*delAzm
     1835                cake = LRAzim[0]+i*delAzm-AzmthOff
    18141836                ind = np.searchsorted(Azm,cake)
    18151837                Plot.plot([arcxI[ind],arcxO[ind]],[arcyI[ind],arcyO[ind]],color='k',dashes=(5,5))
  • trunk/GSASIIstruct.py

    r397 r399  
    109109    return datus[1]
    110110   
    111 def GetHistogramNames(GPXfile):
    112     ''' Returns a list of histogram names found in GSASII gpx file
     111def GetHistogramNames(GPXfile,hType):
     112    """ Returns a list of histogram names found in GSASII gpx file
    113113    input:
    114114        GPXfile = .gpx full file name
     115        hType = list ['PWDR','HKLF']
    115116    return:
    116117        HistogramNames = list of histogram names (types = PWDR & HKLF)
    117     '''
     118    """
    118119    file = open(GPXfile,'rb')
    119120    HistogramNames = []
     
    124125            break
    125126        datum = data[0]
    126         if datum[0][:4] in ['PWDR','HKLF']:
     127        if datum[0][:4] in hType:
    127128            HistogramNames.append(datum[0])
    128129    file.close()
     
    139140    '''
    140141    phaseNames = GetPhaseNames(GPXfile)
     142    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
    141143    phaseData = {}
    142144    for name in phaseNames:
     
    144146    Histograms = {}
    145147    Phases = {}
    146     pId = 0
    147     hId = 0
    148148    for phase in phaseData:
    149149        Phase = phaseData[phase]
    150150        if Phase['Histograms']:
    151151            if phase not in Phases:
     152                pId = phaseNames.index(phase)
    152153                Phase['pId'] = pId
    153                 pId += 1
    154154                Phases[phase] = Phase
    155155            for hist in Phase['Histograms']:
     
    160160                        Histograms[hist] = GetHKLFdata(GPXfile,hist)
    161161                    #future restraint, etc. histograms here           
     162                    hId = histoList.index(hist)
    162163                    Histograms[hist]['hId'] = hId
    163                     hId += 1
    164164    return Histograms,Phases
    165165   
     166def GPXBackup(GPXfile):
     167    import distutils.file_util as dfu
     168    GPXpath,GPXname = ospath.split(GPXfile)
     169    Name = ospath.splitext(GPXname)[0]
     170    files = os.listdir(GPXpath)
     171    last = 0
     172    for name in files:
     173        name = name.split('.')
     174        if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
     175            last = max(last,int(name[1].strip('bak'))+1)
     176    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
     177    dfu.copy_file(GPXfile,GPXback)
     178    return GPXback
     179       
    166180def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData):
    167181    ''' Updates gpxfile from all histograms that are found in any phase
     
    174188    '''
    175189                       
    176     def GPXBackup(GPXfile):
    177         import distutils.file_util as dfu
    178         GPXpath,GPXname = ospath.split(GPXfile)
    179         Name = ospath.splitext(GPXname)[0]
    180         files = os.listdir(GPXpath)
    181         last = 0
    182         for name in files:
    183             name = name.split('.')
    184             if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
    185                 last = max(last,int(name[1].strip('bak'))+1)
    186         GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
    187         dfu.copy_file(GPXfile,GPXback)
    188         return GPXback
    189        
    190190    GPXback = GPXBackup(GPXfile)
    191191    print '\n',135*'-'
     
    200200            break
    201201        datum = data[0]
    202         print 'read: ',datum[0]
     202#        print 'read: ',datum[0]
    203203        if datum[0] == 'Phases':
    204204            for iphase in range(len(data)):
     
    210210        try:
    211211            histogram = Histograms[datum[0]]
    212             print 'found ',datum[0]
     212#            print 'found ',datum[0]
    213213            data[0][1][1] = histogram['Data']
    214214            for datus in data[1:]:
    215                 print '    read: ',datus[0]
     215#                print '    read: ',datus[0]
    216216                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
    217217                    datus[1] = histogram[datus[0]]
     
    223223    outfile.close()
    224224    print 'refinement save successful'
     225   
     226def SetSeqResult(GPXfile,Histograms,SeqResult):
     227    GPXback = GPXBackup(GPXfile)
     228    print '\n',135*'-'
     229    print 'Read from file:',GPXback
     230    print 'Save to file  :',GPXfile
     231    infile = open(GPXback,'rb')
     232    outfile = open(GPXfile,'wb')
     233    while True:
     234        try:
     235            data = cPickle.load(infile)
     236        except EOFError:
     237            break
     238        datum = data[0]
     239        if datum[0] == 'Sequental results':
     240            data[0][1] = SeqResult
     241        try:
     242            histogram = Histograms[datum[0]]
     243            data[0][1][1] = histogram['Data']
     244            for datus in data[1:]:
     245                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
     246                    datus[1] = histogram[datus[0]]
     247        except KeyError:
     248            pass
     249                               
     250        cPickle.dump(data,outfile,1)
     251    infile.close()
     252    outfile.close()
     253    print 'refinement save successful'
     254   
    225255                   
    226256def GetPWDRdata(GPXfile,PWDRname):
     
    878908        pId = Phases[phase]['pId']
    879909        for histogram in HistoPhase:
     910            try:
     911                Histogram = Histograms[histogram]
     912            except KeyError:                       
     913                #skip if histogram not included e.g. in a sequential refinement
     914                continue
    880915            hapData = HistoPhase[histogram]
    881             Histogram = Histograms[histogram]
    882916            hId = Histogram['hId']
    883917            limits = Histogram['Limits'][1]
     
    9681002    return hapVary,hapDict,controlDict
    9691003   
    970 def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms):
     1004def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
    9711005   
    9721006    def PrintSizeAndSig(hapData,sizeSig):
     
    10641098        pId = Phases[phase]['pId']
    10651099        for histogram in HistoPhase:
     1100            try:
     1101                Histogram = Histograms[histogram]
     1102            except KeyError:                       
     1103                #skip if histogram not included e.g. in a sequential refinement
     1104                continue
    10661105            print '\n Phase: ',phase,' in histogram: ',histogram
    10671106            print 130*'-'
    10681107            hapData = HistoPhase[histogram]
    1069             Histogram = Histograms[histogram]
    10701108            hId = Histogram['hId']
    10711109            pfx = str(pId)+':'+str(hId)+':'
     
    10871125                    if pfx+item in sigDict:
    10881126                        PhFrExtPOSig[item] = sigDict[pfx+item]
    1089             if 'Scale' in PhFrExtPOSig:
    1090                 print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
    1091             if 'Extinction' in PhFrExtPOSig:
    1092                 print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
    1093             if hapData['Pref.Ori.'][0] == 'MD':
    1094                 if 'MD' in PhFrExtPOSig:
    1095                     print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
    1096             else:
    1097                 PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
     1127            if Print:
     1128                if 'Scale' in PhFrExtPOSig:
     1129                    print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
     1130                if 'Extinction' in PhFrExtPOSig:
     1131                    print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
     1132                if hapData['Pref.Ori.'][0] == 'MD':
     1133                    if 'MD' in PhFrExtPOSig:
     1134                        print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
     1135                else:
     1136                    PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
    10981137            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
    10991138                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]],
     
    11231162                if pfx+name in sigDict:
    11241163                    SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
    1125             PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
    1126             PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
    1127             PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
     1164            if Print:
     1165                PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
     1166                PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
     1167                PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
    11281168   
    11291169def GetHistogramData(Histograms,Print=True):
     
    12671307    return histVary,histDict,controlDict
    12681308   
    1269 def SetHistogramData(parmDict,sigDict,Histograms):
     1309def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
    12701310   
    12711311    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
     
    13761416            print '\n Histogram: ',histogram,' histogram Id: ',hId
    13771417            print 135*'-'
    1378             print ' Instrument type: ',Sample['Type']
    13791418            print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs'])
    1380             PrintSampleParmsSig(Sample,sampSig)
    1381             PrintInstParmsSig(Inst,instSig)
    1382             PrintBackgroundSig(Background,backSig)
     1419            if Print:
     1420                print ' Instrument type: ',Sample['Type']
     1421                PrintSampleParmsSig(Sample,sampSig)
     1422                PrintInstParmsSig(Inst,instSig)
     1423                PrintBackgroundSig(Background,backSig)
    13831424
    13841425def GetAtomFXU(pfx,FFtables,calcControls,parmDict):
     
    21762217                except IndexError:
    21772218                    pass
    2178     return dMdv   
     2219    return dMdv
     2220
     2221def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
     2222    parmdict.update(zip(varylist,values))
     2223    G2mv.Dict2Map(parmdict)
     2224    Histograms,Phases = HistoPhases
     2225    dMdv = np.empty(0)
     2226    for histogram in Histograms:
     2227        if 'PWDR' in histogram[:4]:
     2228            Histogram = Histograms[histogram]
     2229            hId = Histogram['hId']
     2230            hfx = ':%d:'%(hId)
     2231            Limits = calcControls[hfx+'Limits']
     2232            x,y,w,yc,yb,yd = Histogram['Data']
     2233            xB = np.searchsorted(x,Limits[0])
     2234            xF = np.searchsorted(x,Limits[1])
     2235            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
     2236                varylist,Histogram,Phases,calcControls,pawleyLookup)
     2237            if len(dMdv):
     2238                dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
     2239            else:
     2240                dMdv = dMdvh
     2241    return dMdv
     2242
     2243def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
     2244    parmdict.update(zip(varylist,values))
     2245    Values2Dict(parmdict, varylist, values)
     2246    G2mv.Dict2Map(parmdict)
     2247    Histograms,Phases = HistoPhases
     2248    M = np.empty(0)
     2249    sumwYo = 0
     2250    Nobs = 0
     2251    for histogram in Histograms:
     2252        if 'PWDR' in histogram[:4]:
     2253            Histogram = Histograms[histogram]
     2254            hId = Histogram['hId']
     2255            hfx = ':%d:'%(hId)
     2256            Limits = calcControls[hfx+'Limits']
     2257            x,y,w,yc,yb,yd = Histogram['Data']
     2258            yc *= 0.0                           #zero full calcd profiles
     2259            yb *= 0.0
     2260            yd *= 0.0
     2261            xB = np.searchsorted(x,Limits[0])
     2262            xF = np.searchsorted(x,Limits[1])
     2263            Histogram['Nobs'] = xF-xB
     2264            Nobs += Histogram['Nobs']
     2265            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
     2266            sumwYo += Histogram['sumwYo']
     2267            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
     2268                varylist,Histogram,Phases,calcControls,pawleyLookup)
     2269            yc[xB:xF] += yb[xB:xF]
     2270            yd[xB:xF] = yc[xB:xF]-y[xB:xF]          #yc-yo then all dydv have no '-' needed
     2271            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
     2272            wdy = np.sqrt(w[xB:xF])*(yd[xB:xF])
     2273            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
     2274            M = np.concatenate((M,wdy))
     2275    Histograms['sumwYo'] = sumwYo
     2276    Histograms['Nobs'] = Nobs
     2277    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
     2278    if dlg:
     2279        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile wRp =',Rwp,'%'))[0]
     2280        if not GoOn:
     2281            parmDict['saved values'] = values
     2282            raise Exception         #Abort!!
     2283    return M
     2284   
    21792285                   
    21802286def Refine(GPXfile,dlg):
     
    21822288    import pytexture as ptx
    21832289    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
    2184    
    2185     def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
    2186         parmdict.update(zip(varylist,values))
    2187         G2mv.Dict2Map(parmDict)
    2188         Histograms,Phases = HistoPhases
    2189         dMdv = np.empty(0)
    2190         for histogram in Histograms:
    2191             if 'PWDR' in histogram[:4]:
    2192                 Histogram = Histograms[histogram]
    2193                 hId = Histogram['hId']
    2194                 hfx = ':%d:'%(hId)
    2195                 Limits = calcControls[hfx+'Limits']
    2196                 x,y,w,yc,yb,yd = Histogram['Data']
    2197                 xB = np.searchsorted(x,Limits[0])
    2198                 xF = np.searchsorted(x,Limits[1])
    2199                 dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
    2200                     varylist,Histogram,Phases,calcControls,pawleyLookup)
    2201                 if len(dMdv):
    2202                     dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
    2203                 else:
    2204                     dMdv = dMdvh
    2205         return dMdv
    2206    
    2207     def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
    2208         parmdict.update(zip(varylist,values))
    2209         Values2Dict(parmdict, varylist, values)
    2210         G2mv.Dict2Map(parmDict)
    2211         Histograms,Phases = HistoPhases
    2212         M = np.empty(0)
    2213         sumwYo = 0
    2214         Nobs = 0
    2215         for histogram in Histograms:
    2216             if 'PWDR' in histogram[:4]:
    2217                 Histogram = Histograms[histogram]
    2218                 hId = Histogram['hId']
    2219                 hfx = ':%d:'%(hId)
    2220                 Limits = calcControls[hfx+'Limits']
    2221                 x,y,w,yc,yb,yd = Histogram['Data']
    2222                 yc *= 0.0                           #zero full calcd profiles
    2223                 yb *= 0.0
    2224                 yd *= 0.0
    2225                 xB = np.searchsorted(x,Limits[0])
    2226                 xF = np.searchsorted(x,Limits[1])
    2227                 Histogram['Nobs'] = xF-xB
    2228                 Nobs += Histogram['Nobs']
    2229                 Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
    2230                 sumwYo += Histogram['sumwYo']
    2231                 yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
    2232                     varylist,Histogram,Phases,calcControls,pawleyLookup)
    2233                 yc[xB:xF] += yb[xB:xF]
    2234                 yd[xB:xF] = yc[xB:xF]-y[xB:xF]          #yc-yo then all dydv have no '-' needed
    2235                 Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
    2236                 wdy = np.sqrt(w[xB:xF])*(yd[xB:xF])
    2237                 Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
    2238                 M = np.concatenate((M,wdy))
    2239         Histograms['sumwYo'] = sumwYo
    2240         Histograms['Nobs'] = Nobs
    2241         Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
    2242         if dlg:
    2243             GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile wRp =',Rwp,'%'))[0]
    2244             if not GoOn:
    2245                 parmDict['saved values'] = values
    2246                 raise Exception         #Abort!!
    2247         return M
    22482290   
    22492291    ShowBanner()
     
    23512393#    file.close()
    23522394
     2395def SeqRefine(GPXfile,dlg):
     2396    import cPickle
     2397    import pytexture as ptx
     2398    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
     2399   
     2400    ShowBanner()
     2401    print ' Sequential Refinement'
     2402    G2mv.InitVars()   
     2403    Controls = GetControls(GPXfile)
     2404    ShowControls(Controls)           
     2405    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
     2406    if not Phases:
     2407        print ' *** ERROR - you have no histograms to refine! ***'
     2408        print ' *** Refine aborted ***'
     2409        raise Exception
     2410    if not Histograms:
     2411        print ' *** ERROR - you have no data to refine with! ***'
     2412        print ' *** Refine aborted ***'
     2413        raise Exception
     2414    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables = GetPhaseData(Phases,False)
     2415    if 'Seq Data' in Controls:
     2416        histNames = Controls['Seq Data']
     2417    else:
     2418        histNames = GetHistogramNames(GPXfile,['PWDR',])
     2419    if 'Reverse Seq' in Controls:
     2420        if Controls['Reverse Seq']:
     2421            histNames.reverse()
     2422    SeqResult = {'histNames':histNames}
     2423   
     2424    for ihst,histogram in enumerate(histNames):
     2425        print ihst,histogram
     2426        ifPrint = False
     2427        if dlg:
     2428            dlg.SetTitle('Residual for histogram '+str(ihst))
     2429        calcControls = {}
     2430        calcControls['Natoms'] = Natoms
     2431        calcControls['FFtables'] = FFtables
     2432        varyList = []
     2433        parmDict = {}
     2434        Histo = {histogram:Histograms[histogram],}
     2435        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
     2436        calcControls.update(controlDict)
     2437        histVary,histDict,controlDict = GetHistogramData(Histo,False)
     2438        calcControls.update(controlDict)
     2439        varyList = phaseVary+hapVary+histVary
     2440        if not ihst:
     2441            saveVaryList = varyList[:]
     2442            for i,item in enumerate(saveVaryList):
     2443                items = item.split(':')
     2444                if items[1]:
     2445                    items[1] = ''
     2446                item = ':'.join(items)
     2447                saveVaryList[i] = item
     2448            SeqResult['varyList'] = saveVaryList
     2449        else:
     2450            newVaryList = varyList[:]
     2451            for i,item in enumerate(newVaryList):
     2452                items = item.split(':')
     2453                if items[1]:
     2454                    items[1] = ''
     2455                item = ':'.join(items)
     2456                newVaryList[i] = item
     2457            if newVaryList != SeqResult['varyList']:
     2458                print newVaryList
     2459                print SeqResult['varyList']
     2460                print '**** ERROR - variable list for this histogram does not match previous'
     2461                raise Exception
     2462        parmDict.update(phaseDict)
     2463        parmDict.update(hapDict)
     2464        parmDict.update(histDict)
     2465        GetFprime(calcControls,Histo)
     2466        constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
     2467        groups,parmlist = G2mv.GroupConstraints(constrDict)
     2468        G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
     2469        G2mv.Map2Dict(parmDict,varyList)
     2470   
     2471        while True:
     2472            begin = time.time()
     2473            values =  np.array(Dict2Values(parmDict, varyList))
     2474            Ftol = Controls['min dM/M']
     2475            Factor = Controls['shift factor']
     2476            if Controls['deriv type'] == 'analytic':
     2477                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
     2478                    ftol=Ftol,col_deriv=True,factor=Factor,
     2479                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     2480                ncyc = int(result[2]['nfev']/2)               
     2481            else:           #'numeric'
     2482                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
     2483                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     2484                ncyc = int(result[2]['nfev']/len(varyList))
     2485            runtime = time.time()-begin
     2486            chisq = np.sum(result[2]['fvec']**2)
     2487            Values2Dict(parmDict, varyList, result[0])
     2488            G2mv.Dict2Map(parmDict)
     2489            ApplyXYZshifts(parmDict,varyList)
     2490           
     2491            Rwp = np.sqrt(chisq/Histo['sumwYo'])*100.      #to %
     2492            GOF = chisq/(Histo['Nobs']-len(varyList))
     2493            print '\n Refinement results for histogram: v'+histogram
     2494            print 135*'-'
     2495            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
     2496            print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
     2497            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
     2498            try:
     2499                covMatrix = result[1]*GOF
     2500                sig = np.sqrt(np.diag(covMatrix))
     2501                if np.any(np.isnan(sig)):
     2502                    print '*** Least squares aborted - some invalid esds possible ***'
     2503                    ifPrint = True
     2504                break                   #refinement succeeded - finish up!
     2505            except TypeError:          #result[1] is None on singular matrix
     2506                print '**** Refinement failed - singular matrix ****'
     2507                Ipvt = result[2]['ipvt']
     2508                for i,ipvt in enumerate(Ipvt):
     2509                    if not np.sum(result[2]['fjac'],axis=1)[i]:
     2510                        print 'Removing parameter: ',varyList[ipvt-1]
     2511                        del(varyList[ipvt-1])
     2512                        break
     2513   
     2514        GetFobsSq(Histo,Phases,parmDict,calcControls)
     2515        sigDict = dict(zip(varyList,sig))
     2516        covData = {'variables':result[0],'covMatrix':covMatrix}
     2517        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
     2518        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
     2519        SeqResult[histogram] = covData
     2520    SetSeqResult(GPXfile,Histograms,SeqResult)
     2521
     2522
    23532523def main():
    23542524    arg = sys.argv
Note: See TracChangeset for help on using the changeset viewer.