Changeset 3855


Ignore:
Timestamp:
Mar 15, 2019 7:53:57 PM (4 years ago)
Author:
toby
Message:

major speed up with very big sequential fits

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIIO.py

    r3843 r3855  
    600600        unexpectedObject = True
    601601   
     602def cPickleLoad(fp):
     603    if '2' in platform.python_version_tuple()[0]:
     604        return cPickle.load(fp)
     605    else:
     606       return cPickle.load(fp,encoding='latin-1')
     607
    602608def ProjFileOpen(G2frame,showProvenance=True):
    603609    'Read a GSAS-II project file and load into the G2 data tree'
    … …  
    609615    filep = open(G2frame.GSASprojectfile,'rb')
    610616    if showProvenance: print ('loading from file: '+G2frame.GSASprojectfile)
     617    GPXphase = os.path.splitext(G2frame.GSASprojectfile)[0]+'.seqPhase'
     618    GPXhist = os.path.splitext(G2frame.GSASprojectfile)[0]+'.seqHist'
     619    hist = None
     620    tmpHistIndex = {}
     621    updateFromSeq = False
     622    if os.path.exists(GPXphase) and os.path.exists(GPXhist):
     623        dlg = wx.MessageDialog(G2frame,
     624            'Load results from crashed sequential fit?\nNo deletes the files!', 'Recover partial sequential fit?', wx.YES | wx.NO | wx.CANCEL)
     625        dlg.CenterOnParent()
     626        try:
     627            result = dlg.ShowModal()
     628            if result == wx.ID_OK:
     629                updateFromSeq = True
     630                fp = open(GPXphase,'rb')
     631                data = cPickleLoad(fp) # first block in file should be Phases
     632                if data[0][0] != 'Phases':
     633                    raise Exception('Unexpected block in {} file. How did this happen?'
     634                            .format(GPXphase))
     635                Phases = {}
     636                for name,vals in data[1:]:
     637                    Phases[name] = vals
     638                name,CovData = cPickleLoad(fp)[0] # 2nd block in file should be Covariance
     639                name,RigidBodies = cPickleLoad(fp)[0] # 3rd block in file should be Rigid Bodies
     640                fp.close()
     641                # index the histogram updates
     642                hist = open(GPXhist,'rb')
     643                try:
     644                    while True:
     645                        loc = hist.tell()
     646                        datum = cPickleLoad(hist)[0]
     647                        tmpHistIndex[datum[0]] = loc
     648                except EOFError:
     649                    pass
     650            elif result != wx.ID_CANCEL:
     651                #print('deleting .seqXXXX files')
     652                os.remove(GPXphase)
     653                os.remove(GPXhist)
     654        finally:
     655            dlg.Destroy()
    611656    wx.BeginBusyCursor()
    612657    try:
    613658        while True:
    614659            try:
    615                 if '2' in platform.python_version_tuple()[0]:
    616                     data = cPickle.load(filep)
    617                 else:
    618                     data = cPickle.load(filep,encoding='latin-1')
     660                data = cPickleLoad(filep)
    619661            except EOFError:
    620662                break
    … …  
    629671                #    GSASIIpath.IPyBreak()
    630672            Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=datum[0])
     673            if updateFromSeq and datum[0] == 'Phases':
     674                for pdata in data[1:]:
     675                    if pdata[0] in Phases:
     676                        pdata[1].update(Phases[pdata[0]])
     677            elif updateFromSeq and datum[0] == 'Covariance':
     678                data[0][1] = CovData
     679            elif updateFromSeq and datum[0] == 'Rigid bodies':
     680                data[0][1] = RigidBodies
     681            elif updateFromSeq and datum[0] in tmpHistIndex:
     682                hist.seek(tmpHistIndex[datum[0]])
     683                hdata = cPickleLoad(hist)
     684                xferItems = ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']
     685                hItems = {name:j+1 for j,(name,val) in enumerate(hdata[1:]) if name in xferItems}
     686                for j,(name,val) in enumerate(data[1:]):
     687                    if name not in xferItems: continue
     688                    data[j+1][1] = hdata[hItems[name]][1]
    631689            if datum[0].startswith('PWDR'):               
    632690                if 'ranId' not in datum[1][0]: # patch: add random Id if not present
    … …  
    693751        wx.EndBusyCursor()
    694752        G2frame.Status.SetStatusText('Mouse RB drag/drop to reorder',0)
     753    if hist:
     754        hist.close()
     755        #print('deleting .seqXXXX files')
     756        os.remove(GPXphase)
     757        os.remove(GPXhist)
    695758    G2frame.SetTitleByGPX()
    696759   

  • trunk/GSASIIdataGUI.py

    r3848 r3855  
    38163816            else:
    38173817                pth = '.'
     3818            #if GSASIIpath.GetConfigValue('debug'): print('debug: open from '+pth)
    38183819            dlg = wx.FileDialog(self, 'Choose GSAS-II project file', pth,
    38193820                wildcard='GSAS-II project file (*.gpx)|*.gpx',style=wx.FD_OPEN)
    … …  
    38343835        if not filename:
    38353836            GetGPX()
     3837            filename = self.GSASprojectfile
    38363838        else:
    38373839            try:
    … …  
    38783880            self.GPXtree.Expand(Id)
    38793881            SelectDataTreeItem(self,Id)
    3880 #            self.GPXtree.SelectItem(Id)
    38813882        elif phaseId:
    3882 #            self.GPXtree.SelectItem(phaseId)
    38833883            SelectDataTreeItem(self,phaseId)
    38843884        self.CheckNotebook()
    … …  
    43604360                if name[:4] in hType:
    43614361                    HistogramNames.append(name)       
    4362                 item, cookie = self.GPXtree.GetNextChild(self.root, cookie)               
    4363 
     4362                item, cookie = self.GPXtree.GetNextChild(self.root, cookie)
    43644363        return HistogramNames
     4364   
     4365    def GetHistogramNamesID(self,hType):
     4366        """ Returns a list of histogram names found in the GSASII data tree
     4367        and a lookup table of their Id values. Should replace GetHistogramNames
     4368        since that will not be much faster (and there may be real speed gains from
     4369        caching the Ids rather than keep searching for them).
     4370
     4371        N.B routine :func:`GSASIIstrIO.GetHistogramNames` also exists to
     4372        get same info, but from GPX file.
     4373       
     4374        :param str hType: list of histogram types
     4375        :return: list of histogram names and a dict of histogram Ids
     4376           keyed by histogram name.
     4377        """
     4378        HistogramNames = []
     4379        HistogramIds = {}
     4380        if self.GPXtree.GetCount():
     4381            item, cookie = self.GPXtree.GetFirstChild(self.root)
     4382            while item:
     4383                name = self.GPXtree.GetItemText(item)
     4384                if name[:4] in hType:
     4385                    HistogramNames.append(name)       
     4386                    HistogramIds[name] = item
     4387                item, cookie = self.GPXtree.GetNextChild(self.root, cookie)
     4388        return HistogramNames,HistogramIds
    43654389                   
    43664390    def GetUsedHistogramsAndPhasesfromTree(self):
    … …  
    43814405        phaseNames = self.GetPhaseNames()
    43824406        phaseData = self.GetPhaseData()
    4383         histoList = self.GetHistogramNames(['PWDR','HKLF'])
     4407        histoList,histIdList = self.GetHistogramNamesID(['PWDR','HKLF'])
    43844408
    43854409        for phase in phaseData:
    … …  
    43944418                        Phases[phase] = Phase
    43954419                    if hist not in Histograms and Phase['Histograms'][hist]['Use']:
    4396                         item = GetGPXtreeItemId(self,self.root,hist)
     4420                        item = histIdList[hist]
    43974421                        if item:
    43984422                            if 'PWDR' in hist[:4]:
    … …  
    45134537        Called from the Calculate/Refine menu.
    45144538        '''
    4515         G2cnstG.CheckAllScalePhaseFractions(self)
    45164539        if self.testSeqRefineMode():
    45174540            self.OnSeqRefine(event)
    45184541            return
    4519         # Id = GetGPXtreeItemId(self,self.root,'Sequential results')
    4520         # if Id:
    4521         #     dlg = wx.MessageDialog(
    4522         #         self,
    4523         #         'Your last refinement was sequential. Continue with "Refine", removing previous sequential results?',
    4524         #         'Remove sequential results?',wx.OK|wx.CANCEL)
    4525         #     if dlg.ShowModal() == wx.ID_OK:
    4526         #         self.GPXtree.Delete(Id)
    4527         #         dlg.Destroy()
    4528         #     else:
    4529         #         dlg.Destroy()
    4530         #         return
     4542        G2cnstG.CheckAllScalePhaseFractions(self) # can be slow for sequential fits, skip
    45314543        self.OnFileSave(event)
    45324544        # check that constraints are OK here
    … …  
    46404652    def OnSeqRefine(self,event):
    46414653        '''Perform a sequential refinement.
    4642         Called from self.OnRefine (Called from the Calculate/Refine menu)
    4643        
    4644         temporarily called from the Calculate/Sequential refine menu (to be removed)
     4654        Called from self.OnRefine (Which is called from the Calculate/Refine menu)
    46454655        '''
    46464656        seqList = self.testSeqRefineMode()

  • trunk/GSASIIpwdGUI.py

    r3848 r3855  
    44904490        if phaseName not in ['Unknown',]:
    44914491            pId = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Phases')
    4492             phaseId =  G2gd.GetGPXtreeItemId(G2frame,pId,phaseName)
    4493             if not phaseId:         #phase deleted
    4494                 return None
    4495             General = G2frame.GPXtree.GetItemPyData(phaseId)['General']
    4496             if General.get('Modulated',False):
    4497                 Super = 1
     4492            if pId: # phase section missing from file (unusual)
     4493                phaseId =  G2gd.GetGPXtreeItemId(G2frame,pId,phaseName)
     4494                if phaseId:         #is phase deleted?
     4495                    General = G2frame.GPXtree.GetItemPyData(phaseId)['General']
     4496                    if General.get('Modulated',False):
     4497                        Super = 1
    44984498        rowLabels = []
    44994499        if HKLF:

  • trunk/GSASIIstrIO.py

    r3849 r3855  
    4949ateln2 = 8.0*math.log(2.0)
    5050
     51#===============================================================================
     52# Support for GPX file reading
     53#===============================================================================
    5154def cPickleLoad(fp):
    5255    if '2' in platform.python_version_tuple()[0]:
    … …  
    5558       return cPickle.load(fp,encoding='latin-1')
    5659
     60gpxIndex = {}; gpxNamelist = []; gpxSize = -1
     61'''Global variables used in :func:`IndexGPX` to see if file has changed
     62(gpxSize) and to index where to find each 1st-level tree item in the file.
     63'''
     64tmpHistIndex = {}
     65'Global variable used to index contents of .seqHist file'
     66
    5767def GetFullGPX(GPXfile):
    5868    ''' Returns complete contents of GSASII gpx file.
    … …  
    6777      * nameList (list) has names of main tree entries & subentries used to reconstruct project file
    6878    '''
     79    return IndexGPX(GPXfile,read=True)
     80
     81def IndexGPX(GPXfile,read=False):
     82    '''Create an index to a GPX file, optionally the file into memory.
     83    The byte size of the GPX file is saved. If this routine is called
     84    again, and if this size does not change, indexing is not repeated
     85    since it is assumed the file has not changed (this can be overriden
     86    by setting read=True).
     87
     88    :param str GPXfile: full .gpx file name
     89    :returns: Project,nameList if read=, where
     90
     91      * Project (dict) is a representation of gpx file following the GSAS-II
     92        tree structure for each item: key = tree name (e.g. 'Controls',
     93        'Restraints', etc.), data is dict
     94      * nameList (list) has names of main tree entries & subentries used to reconstruct project file
     95    '''
     96    global gpxSize
     97    if gpxSize == os.path.getsize(GPXfile) and not read:
     98        return
     99    global gpxIndex
     100    gpxIndex = {}
     101    global gpxNamelist
     102    gpxNamelist = []
     103    if GSASIIpath.GetConfigValue('debug'): print("DBG: Indexing GPX file")
     104    gpxSize = os.path.getsize(GPXfile)
    69105    fp = open(GPXfile,'rb')
    70106    Project = {}
    71     nameList = []
    72107    try:
    73108        while True:
    74             try:
    75                 data = cPickleLoad(fp)
    76             except EOFError:
    77                 break
     109            pos = fp.tell()
     110            data = cPickleLoad(fp)
    78111            datum = data[0]
    79             Project[datum[0]] = {'data':datum[1]}
    80             nameList.append([datum[0],])
     112            gpxIndex[datum[0]] = pos
     113            if read: Project[datum[0]] = {'data':datum[1]}
     114            gpxNamelist.append([datum[0],])
    81115            for datus in data[1:]:
    82                 Project[datum[0]][datus[0]] = datus[1]
    83                 nameList[-1].append(datus[0])
     116                if read: Project[datum[0]][datus[0]] = datus[1]
     117                gpxNamelist[-1].append(datus[0])
    84118        # print('project load successful')
     119    except EOFError:
     120        pass
    85121    except Exception as msg:
    86122        print('Read Error:',msg)
    … …  
    88124    finally:
    89125        fp.close()
    90     return Project,nameList
     126    if read: return Project,gpxNamelist   
    91127   
    92128def GetControls(GPXfile):
    … …  
    97133    '''
    98134    Controls = copy.copy(G2obj.DefaultControls)
    99     fl = open(GPXfile,'rb')
    100     while True:
    101         try:
    102             data = cPickleLoad(fl)
    103         except EOFError:
    104             break
    105         datum = data[0]
    106         if datum[0] == 'Controls':
    107             Controls.update(datum[1])
    108     fl.close()
     135    IndexGPX(GPXfile)
     136    pos = gpxIndex.get('Controls')
     137    if pos is None:
     138        print('Warning: Controls not found in gpx file {}'.format(GPXfile))
     139        return Controls
     140    fp = open(GPXfile,'rb')
     141    fp.seek(pos)
     142    datum = cPickleLoad(fp)[0]
     143    fp.close()
     144    Controls.update(datum[1])
    109145    return Controls
    110    
     146
    111147def GetConstraints(GPXfile):
    112148    '''Read the constraints from the GPX file and interpret them
    … …  
    115151    and :func:`GSASIIstrMain.SeqRefine`.
    116152    '''
     153    IndexGPX(GPXfile)
    117154    fl = open(GPXfile,'rb')
    118     while True:
    119         try:
    120             data = cPickleLoad(fl)
    121         except EOFError:
    122             break
    123         datum = data[0]
    124         if datum[0] == 'Constraints':
    125             constList = []
    126             for item in datum[1]:
    127                 if item.startswith('_'): continue
    128                 constList += datum[1][item]
    129             fl.close()
    130             constDict,fixedList,ignored = ProcessConstraints(constList)
    131             if ignored:
    132                 print ('%d Constraints were rejected. Was a constrained phase, histogram or atom deleted?'%ignored)
    133             return constDict,fixedList
     155    pos = gpxIndex.get('Constraints')
     156    if pos is None:
     157        raise Exception("No constraints in GPX file")
     158    fl.seek(pos)
     159    datum = cPickleLoad(fl)[0]
    134160    fl.close()
    135     raise Exception("No constraints in GPX file")
     161    constList = []
     162    for item in datum[1]:
     163        if item.startswith('_'): continue
     164        constList += datum[1][item]
     165    constDict,fixedList,ignored = ProcessConstraints(constList)
     166    if ignored:
     167        print ('{} Constraints were rejected. Was a constrained phase, histogram or atom deleted?'.format(ignored))
     168    return constDict,fixedList
    136169   
    137170def ProcessConstraints(constList):
    … …  
    226259    '''Load constraints and related info and return any error or warning messages
    227260    This is done from the GPX file rather than the tree.
    228     This is called before a refinement is launched (OnRefine and OnSeqRefine),
    229     where the tree could be used.
    230    
     261
    231262    :param dict seqHist: defines a specific histogram to be loaded for a sequential
    232263       refinement, if None (default) all are loaded.
    … …  
    247278    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
    248279        GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False) # generates atom symmetry constraints
    249     hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False,resetRefList=True)
     280    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False,resetRefList=False)
    250281    histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
    251282    varyList = rbVary+phaseVary+hapVary+histVary
    … …  
    282313    Throws an exception if not found in the .GPX file
    283314    '''
     315    IndexGPX(GPXfile)
    284316    fl = open(GPXfile,'rb')
    285     while True:
    286         try:
    287             data = cPickleLoad(fl)
    288         except EOFError:
    289             break
    290         datum = data[0]
    291         if datum[0] == 'Restraints':
    292             restraintDict = datum[1]
     317    pos = gpxIndex.get('Restraints')
     318    if pos is None:
     319        raise Exception("No Restraints in GPX file")
     320    fl.seek(pos)
     321    datum = cPickleLoad(fl)[0]
    293322    fl.close()
    294     return restraintDict
     323    return datum[1]
    295324   
    296325def GetRigidBodies(GPXfile):
    297326    '''Read the rigid body models from the GPX file
    298327    '''
     328    IndexGPX(GPXfile)
    299329    fl = open(GPXfile,'rb')
    300     while True:
    301         try:
    302             data = cPickleLoad(fl)
    303         except EOFError:
    304             break
    305         datum = data[0]
    306         if datum[0] == 'Rigid bodies':
    307             rigidbodyDict = datum[1]
     330    pos = gpxIndex.get('Rigid bodies')
     331    if pos is None:
     332        raise Exception("No Rigid bodies in GPX file")
     333    fl.seek(pos)
     334    datum = cPickleLoad(fl)[0]
    308335    fl.close()
    309     return rigidbodyDict
     336    return datum[1]
    310337       
    311338def GetFprime(controlDict,Histograms):
    … …  
    348375    :return: list of phase names
    349376    '''
     377    IndexGPX(GPXfile)
    350378    fl = open(GPXfile,'rb')
    351     PhaseNames = []
    352     while True:
    353         try:
    354             data = cPickleLoad(fl)
    355         except EOFError:
    356             break
    357         datum = data[0]
    358         if 'Phases' == datum[0]:
    359             for datus in data[1:]:
    360                 PhaseNames.append(datus[0])
     379    pos = gpxIndex.get('Phases')
     380    if pos is None:
     381        raise Exception("No Phases in GPX file")
     382    fl.seek(pos)
     383    data = cPickleLoad(fl)
    361384    fl.close()
    362     return PhaseNames
     385    return [datus[0] for datus in data[1:]]
    363386
    364387def GetAllPhaseData(GPXfile,PhaseName):
    … …  
    367390    :param str GPXfile: full .gpx file name
    368391    :param str PhaseName: phase name
    369     :return: phase dictionary
     392    :return: phase dictionary or None if PhaseName is not present
    370393    '''       
     394    IndexGPX(GPXfile)
    371395    fl = open(GPXfile,'rb')
    372     while True:
    373         try:
    374             data = cPickleLoad(fl)
    375         except EOFError:
    376             break
    377         datum = data[0]
    378         if 'Phases' == datum[0]:
    379             for datus in data[1:]:
    380                 if datus[0] == PhaseName:
    381                     break
     396    pos = gpxIndex.get('Phases')
     397    if pos is None:
     398        raise Exception("No Phases in GPX file")
     399    fl.seek(pos)
     400    data = cPickleLoad(fl)
    382401    fl.close()
    383     return datus[1]
     402
     403    for datus in data[1:]:
     404        if datus[0] == PhaseName:
     405            return datus[1]
    384406   
    385407def GetHistograms(GPXfile,hNames):
    … …  
    391413
    392414    """
     415    IndexGPX(GPXfile)
    393416    fl = open(GPXfile,'rb')
    394417    Histograms = {}
    395     while True:
    396         try:
    397             data = cPickleLoad(fl)
    398         except EOFError:
    399             break
     418    for hist in hNames:
     419        pos = gpxIndex.get(hist)
     420        if pos is None:
     421            raise Exception("Histogram {} not found in GPX file".format(hist))
     422        fl.seek(pos)
     423        data = cPickleLoad(fl)
    400424        datum = data[0]
    401         hist = datum[0]
    402         if hist in hNames:
    403             if 'PWDR' in hist[:4]:
    404                 PWDRdata = {}
    405                 PWDRdata.update(datum[1][0])        #weight factor
    406                 PWDRdata['Data'] = ma.array(ma.getdata(datum[1][1]))          #masked powder data arrays/clear previous masks
    407                 PWDRdata[data[2][0]] = data[2][1]       #Limits & excluded regions (if any)
    408                 PWDRdata[data[3][0]] = data[3][1]       #Background
    409                 PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
    410                 PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
    411                 try:
    412                     PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
    413                 except IndexError:
    414                     PWDRdata['Reflection Lists'] = {}
    415                 PWDRdata['Residuals'] = {}
    416    
    417                 Histograms[hist] = PWDRdata
    418             elif 'HKLF' in hist[:4]:
    419                 HKLFdata = {}
    420                 HKLFdata.update(datum[1][0])        #weight factor
     425        if 'PWDR' in hist[:4]:
     426            PWDRdata = {}
     427            PWDRdata.update(datum[1][0])        #weight factor
     428            PWDRdata['Data'] = ma.array(ma.getdata(datum[1][1]))          #masked powder data arrays/clear previous masks
     429            PWDRdata[data[2][0]] = data[2][1]       #Limits & excluded regions (if any)
     430            PWDRdata[data[3][0]] = data[3][1]       #Background
     431            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
     432            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
     433            try:
     434               PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
     435            except IndexError:
     436                PWDRdata['Reflection Lists'] = {}
     437            PWDRdata['Residuals'] = {}
     438            Histograms[hist] = PWDRdata
     439        elif 'HKLF' in hist[:4]:
     440            HKLFdata = {}
     441            HKLFdata.update(datum[1][0])        #weight factor
    421442#patch
    422                 if 'list' in str(type(datum[1][1])):
    423                 #if isinstance(datum[1][1],list):
    424                     RefData = {'RefList':[],'FF':{}}
    425                     for ref in datum[1][1]:
    426                         RefData['RefList'].append(ref[:11]+[ref[13],])
    427                     RefData['RefList'] = np.array(RefData['RefList'])
    428                     datum[1][1] = RefData
     443            if 'list' in str(type(datum[1][1])):
     444            #if isinstance(datum[1][1],list):
     445                RefData = {'RefList':[],'FF':{}}
     446                for ref in datum[1][1]:
     447                    RefData['RefList'].append(ref[:11]+[ref[13],])
     448                RefData['RefList'] = np.array(RefData['RefList'])
     449                datum[1][1] = RefData
    429450#end patch
    430                 datum[1][1]['FF'] = {}
    431                 HKLFdata['Data'] = datum[1][1]
    432                 HKLFdata[data[1][0]] = data[1][1]       #Instrument parameters
    433                 HKLFdata['Reflection Lists'] = None
    434                 HKLFdata['Residuals'] = {}
    435                 Histograms[hist] = HKLFdata           
     451            datum[1][1]['FF'] = {}
     452            HKLFdata['Data'] = datum[1][1]
     453            HKLFdata[data[1][0]] = data[1][1]       #Instrument parameters
     454            HKLFdata['Reflection Lists'] = None
     455            HKLFdata['Residuals'] = {}
     456            Histograms[hist] = HKLFdata           
    436457    fl.close()
    437458    return Histograms
    438459   
    439 def GetHistogramNames(GPXfile,hType):
    440     """ Returns a list of histogram names found in GSASII gpx file
     460def GetHistogramNames(GPXfile,hTypes):
     461    """ Returns a list of histogram names found in a GSAS-II .gpx file that
     462    match specifed histogram types. Names are returned in the order they
     463    appear in the file.
    441464
    442465    :param str GPXfile: full .gpx file name
    443     :param str hType: list of histogram types
     466    :param str hTypes: list of histogram types
    444467    :return: list of histogram names (types = PWDR & HKLF)
    445468
    446469    """
    447     fl = open(GPXfile,'rb')
    448     HistogramNames = []
    449     while True:
    450         try:
    451             data = cPickleLoad(fl)
    452         except EOFError:
    453             break
    454         datum = data[0]
    455         if datum[0][:4] in hType:
    456             HistogramNames.append(datum[0])
    457     fl.close()
    458     return HistogramNames
     470    IndexGPX(GPXfile)
     471    return [n[0] for n in gpxNamelist if n[0][:4] in hTypes]
    459472   
    460473def GetUsedHistogramsAndPhases(GPXfile):
    … …  
    529542def GPXBackup(GPXfile,makeBack=True):
    530543    '''
    531     makes a backup of the current .gpx file (?)
     544    makes a backup of the specified .gpx file
    532545   
    533546    :param str GPXfile: full .gpx file name
    … …  
    553566    ''' Updates gpxfile from all histograms that are found in any phase
    554567    and any phase that used a histogram. Also updates rigid body definitions.
    555 
     568    This is used for non-sequential fits, but not for sequential fitting.
    556569
    557570    :param str GPXfile: full .gpx file name
    … …  
    617630    :returns: a dict containing the sequential results table
    618631    '''
     632    IndexGPX(GPXfile)
     633    pos = gpxIndex.get('Sequential results')
     634    if pos is None:
     635        return {}
    619636    fl = open(GPXfile,'rb')
    620     SeqResult = {}
    621     while True:
    622         try:
    623             data = cPickleLoad(fl)
    624         except EOFError:
     637    fl.seek(pos)
     638    datum = cPickleLoad(fl)[0]
     639    fl.close()
     640    return datum[1]
     641   
     642def SetupSeqSavePhases(GPXfile):
     643    '''Initialize the files used to save intermediate results from
     644    sequential fits.
     645    '''
     646    IndexGPX(GPXfile)
     647    # load initial Phase results from GPX
     648    fl = open(GPXfile,'rb')
     649    pos = gpxIndex.get('Phases')
     650    if pos is None:
     651        raise Exception("No Phases in GPX file")
     652    fl.seek(pos)
     653    data = cPickleLoad(fl)
     654    fl.close()
     655    # create GPX-like file to store latest Phase info; init with start vals
     656    GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase'
     657    fp = open(GPXphase,'wb')
     658    cPickle.dump(data,fp,1)
     659    fp.close()
     660    # create empty file for histogram info
     661    tmpHistIndex.clear()
     662    GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist'
     663    fp = open(GPXhist,'wb')
     664    fp.close()
     665
     666def SaveUpdatedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData):
     667    '''
     668    Save phase and histogram information into "pseudo-gpx" files. The phase
     669    information is overwritten each time this is called, but histogram information is
     670    appended after each sequential step.
     671
     672    :param str GPXfile: full .gpx file name
     673    :param dict Histograms: dictionary of histograms as {name:data,...}
     674    :param dict Phases: dictionary of phases that use histograms
     675    :param dict RigidBodies: dictionary of rigid bodies
     676    :param dict CovData: dictionary of refined variables, varyList, & covariance matrix
     677    '''
     678                       
     679    import distutils.file_util as dfu
     680   
     681    GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase'
     682    fp = open(GPXphase,'rb')
     683    data = cPickleLoad(fp) # first block in file should be Phases
     684    if data[0][0] != 'Phases':
     685        raise Exception('Unexpected block in {} file. How did this happen?'
     686                            .format(GPXphase))
     687    fp.close()
     688    # update previous phase info
     689    for datum in data[1:]:
     690        if datum[0] in Phases:
     691            datum[1].update(Phases[datum[0]])
     692    # save latest Phase/refinement info
     693    fp = open(GPXphase,'wb')
     694    cPickle.dump(data,fp,1)
     695    cPickle.dump([['Covariance',CovData]],fp,1)
     696    cPickle.dump([['Rigid bodies',RigidBodies]],fp,1)
     697    fp.close()
     698    # create an entry that looks like a PWDR tree item
     699    for key in Histograms:
     700        if key.startswith('PWDR '):
    625701            break
    626         datum = data[0]
    627         if datum[0] == 'Sequential results':
    628             SeqResult = datum[1]
    629     fl.close()
    630     return SeqResult
     702    else:
     703        raise Exception('No PWDR entry in Histogram dict!')
     704    histname = key
     705    hist = copy.deepcopy(Histograms[key])
     706    xfer_dict = {'Index Peak List': [[], []],
     707                   'Comments': [],
     708                   'Unit Cells List': [],
     709                   'Peak List': {'peaks': [], 'sigDict': {}},
     710                   }
     711    histData = hist['Data']
     712    del hist['Data']
     713    for key in ('Limits','Background','Instrument Parameters',
     714                    'Sample Parameters','Reflection Lists'):
     715        xfer_dict[key] = hist[key]
     716        del hist[key]
     717    # xform into a gpx-type entry
     718    data = []
     719    data.append([histname,[hist,histData,histname]])       
     720    for key in ['Comments','Limits','Background','Instrument Parameters',
     721             'Sample Parameters','Peak List','Index Peak List',
     722             'Unit Cells List','Reflection Lists']:
     723        data.append([key,xfer_dict[key]])
     724    # append histogram to histogram info
     725    GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist'
     726    fp = open(GPXhist,'ab')
     727    tmpHistIndex[histname] = fp.tell()
     728    cPickle.dump(data,fp,1)
     729    fp.close()
     730    return   
    631731   
    632732def SetSeqResult(GPXfile,Histograms,SeqResult):
    … …  
    641741    print ('Read from file:'+GPXback)
    642742    print ('Save to file  :'+GPXfile)
     743    GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase'
     744    fp = open(GPXphase,'rb')
     745    data = cPickleLoad(fp) # first block in file should be Phases
     746    if data[0][0] != 'Phases':
     747        raise Exception('Unexpected block in {} file. How did this happen?'
     748                            .format(GPXphase))
     749    Phases = {}
     750    for name,vals in data[1:]:
     751        Phases[name] = vals
     752    name,CovData = cPickleLoad(fp)[0] # 2nd block in file should be Covariance
     753    name,RigidBodies = cPickleLoad(fp)[0] # 3rd block in file should be Rigid Bodies
     754    fp.close()
     755    GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist'
     756    hist = open(GPXhist,'rb')
    643757    infile = open(GPXback,'rb')
    644758    outfile = open(GPXfile,'wb')
    … …  
    651765        if datum[0] == 'Sequential results':
    652766            data[0][1] = SeqResult
    653         # reset the Copy Next flag, since it should not be needed twice in a row
    654         if datum[0] == 'Controls':
     767        elif datum[0] == 'Phases':
     768            for pdata in data[1:]:
     769                if pdata[0] in Phases:
     770                    pdata[1].update(Phases[pdata[0]])
     771        elif datum[0] == 'Covariance':
     772            data[0][1] = CovData
     773        elif datum[0] == 'Rigid bodies':
     774            data[0][1] = RigidBodies
     775        elif datum[0] == 'Controls': # reset the Copy Next flag after a sequential fit
    655776            data[0][1]['Copy2Next'] = False
    656         try:
    657             histogram = Histograms[datum[0]]
    658             data[0][1][1] = list(histogram['Data'])
    659             for datus in data[1:]:
    660                 if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
    661                     datus[1] = histogram[datus[0]]
    662         except KeyError:
    663             pass
    664                                
     777        elif datum[0] in tmpHistIndex:
     778            hist.seek(tmpHistIndex[datum[0]])
     779            hdata = cPickleLoad(hist)
     780            xferItems = ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']
     781            hItems = {name:j+1 for j,(name,val) in enumerate(hdata[1:]) if name in xferItems}
     782            for j,(name,val) in enumerate(data[1:]):
     783                if name not in xferItems: continue
     784                data[j+1][1] = hdata[hItems[name]][1]
     785            # old code update from Histograms array
     786            #histogram = Histograms[datum[0]]
     787            #data[0][1][1] = list(histogram['Data'])
     788            #for datus in data[1:]:
     789            #    if datus[0] in
     790            #        datus[1] = histogram[datus[0]]
    665791        cPickle.dump(data,outfile,1)
    666792    infile.close()
    667793    outfile.close()
    668     print ('GPX file save successful')
    669                        
     794    # clean up tmp files
     795    tmpHistIndex.clear()
     796    os.remove(GPXphase)
     797    os.remove(GPXhist)
     798    print ('GPX file merge completed')
     799
     800#==============================================================================
     801# Refinement routines
     802#==============================================================================
    670803def ShowBanner(pFile=None):
    671804    'Print authorship, copyright and citation notice'

  • trunk/GSASIIstrMain.py

    r3845 r3855  
    364364    Histo = {}
    365365    NewparmDict = {}
     366    G2stIO.SetupSeqSavePhases(GPXfile)
    366367    for ihst,histogram in enumerate(histNames):
    367368        if GSASIIpath.GetConfigValue('Show_timing'): t1 = time.time()
    … …  
    415416        # do constraint processing
    416417        #reload(G2mv) # debug
    417         if GSASIIpath.GetConfigValue('Show_timing'):
    418             t2 = time.time()
    419             print("#1 show time {:.2f} sec.".format(t2-t1))
    420             t1 = t2
    421418        G2mv.InitVars()
    422419        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
    423         if GSASIIpath.GetConfigValue('Show_timing'):
    424             t2 = time.time()
    425             print("#2 show time {:.2f} sec.".format(t2-t1))
    426             t1 = t2
    427420        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
    428421        msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
    … …  
    484477        printFile.write('\n Refinement results for histogram: %s\n'%histogram)
    485478        printFile.write(135*'-'+'\n')
    486         if GSASIIpath.GetConfigValue('Show_timing'):
    487             t2 = time.time()
    488             print("#3 show time {:.2f} sec.".format(t2-t1))
    489             t1 = t2
    490479        if True:
    491480#        try:
    492481            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
    493482                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg)
    494             if GSASIIpath.GetConfigValue('Show_timing'):
    495                 t2 = time.time()
    496                 print("#4a show time {:.2f} sec.".format(t2-t1))
    497                 t1 = t2
    498483            if PlotFunction:
    499484                PlotFunction(G2frame,Histo[histogram]['Data'],histogram)
    500 
    501             if GSASIIpath.GetConfigValue('Show_timing'):
    502                 t2 = time.time()
    503                 print("#4b show time {:.2f} sec.".format(t2-t1))
    504                 t1 = t2
    505485            print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
    506486                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2']))
    … …  
    525505                'parmDict':parmDict}
    526506            SeqResult[histogram] = histRefData
    527             if GSASIIpath.GetConfigValue('Show_timing'):
    528                 t2 = time.time()
    529                 print("#4c show time {:.2f} sec.".format(t2-t1))
    530                 t1 = t2
    531507            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
    532             if GSASIIpath.GetConfigValue('Show_timing'):
    533                 t2 = time.time()
    534                 print("#5a show time {:.2f} sec.".format(t2-t1))
    535                 t1 = t2
    536     #        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
     508#            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
    537509            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
    538             if GSASIIpath.GetConfigValue('Show_timing'):
    539                 t2 = time.time()
    540                 print("#5b show time {:.2f} sec.".format(t2-t1))
    541                 t1 = t2
    542510            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile)
    543             if GSASIIpath.GetConfigValue('Show_timing'):
    544                 t2 = time.time()
    545                 print("#5c show time {:.2f} sec.".format(t2-t1))
    546                 t1 = t2
    547             G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack)
    548             if GSASIIpath.GetConfigValue('Show_timing'):
    549                 t2 = time.time()
    550                 print("#5d show time {:.2f} sec.".format(t2-t1))
    551                 t1 = t2
     511            G2stIO.SaveUpdatedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData)
    552512            makeBack = False
    553513            NewparmDict = {}
    … …  
    570530        if GSASIIpath.GetConfigValue('Show_timing'):
    571531            t2 = time.time()
    572             print("#6 show time {:.2f} sec.".format(t2-t1))
     532            print("Fit step time {:.2f} sec.".format(t2-t1))
    573533            t1 = t2
    574534    SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
Note: See TracChangeset for help on using the changeset viewer.