source: trunk/GSASIIscriptable.py @ 3023

Last change on this file since 3023 was 3023, checked in by odonnell, 4 years ago

remove wx dependency for GSASIIscriptable

File size: 82.3 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3########### SVN repository information ###################
4# $Date: 2017-04-12 15:12:45 -0500 (Wed, 12 Apr 2017) $
5# $Author: vondreele $
6# $Revision: 2777 $
7# $URL: https://subversion.xray.aps.anl.gov/pyGSAS/trunk/GSASIIscriptable.py $
8# $Id: GSASIIIO.py 2777 2017-04-12 20:12:45Z vondreele $
9########### SVN repository information ###################
10"""
11*GSASIIscriptable: Scripting Tools*
12-----------------------------------
13
14Routines for reading, writing, modifying and creating GSAS-II project (.gpx) files.
15
16Supports a command line interface as well.
17
18Look at :class:`G2Project` to start.
19
20=====================
21Refinement parameters
22=====================
23
24There are three classes of refinement parameters:
25
26    * Histogram. Turned on and off through :func:`~G2PwdrData.set_refinements`
27      and :func:`~G2PwdrData.clear_refinements`
28    * Phase. Turned on and off through :func:`~G2Phase.set_refinements`
29      and :func:`~G2Phase.clear_refinements`
30    * Histogram-and-phase (HAP). Turned on and off through
31      :func:`~G2Phase.set_HAP_refinements` and :func:`~G2Phase.clear_HAP_refinements`
32"""
33from __future__ import division, print_function # needed?
34import argparse
35import os.path as ospath
36import datetime as dt
37import sys
38import cPickle
39import imp
40import copy
41import os
42import random as ran
43
44import numpy.ma as ma
45import scipy.interpolate as si
46import numpy as np
47import scipy as sp
48
49import GSASIIpath
50GSASIIpath.SetBinaryPath(False) # would rather have this in __name__ == '__main__' stanza
51import GSASIIobj as G2obj
52import GSASIIpwd as G2pwd
53import GSASIIstrMain as G2strMain
54import GSASIIspc as G2spc
55import GSASIIElem as G2elem
56
57# Delay imports to not slow down small scripts
58G2fil = None
59PwdrDataReaders = []
60PhaseReaders = []
61
62
63def LoadG2fil():
64    """Delay importing this module, it is slow"""
65    global G2fil
66    global PwdrDataReaders
67    global PhaseReaders
68    if G2fil is None:
69        import GSASIIfiles
70        G2fil = GSASIIfiles
71        PwdrDataReaders = G2fil.LoadImportRoutines("pwd", "Powder_Data")
72        PhaseReaders = G2fil.LoadImportRoutines("phase", "Phase")
73
74
75def LoadDictFromProjFile(ProjFile):
76    '''Read a GSAS-II project file and load items to dictionary
77    :param str ProjFile: GSAS-II project (name.gpx) full file name
78    :returns: Project,nameList, where
79
80      * Project (dict) is a representation of gpx file following the GSAS-II tree struture
81        for each item: key = tree name (e.g. 'Controls','Restraints',etc.), data is dict
82        data dict = {'data':item data whch may be list, dict or None,'subitems':subdata (if any)}
83      * nameList (list) has names of main tree entries & subentries used to reconstruct project file
84
85    Example for fap.gpx::
86
87      Project = {                 #NB:dict order is not tree order
88        u'Phases':{'data':None,'fap':{phase dict}},
89        u'PWDR FAP.XRA Bank 1':{'data':[histogram data list],'Comments':comments,'Limits':limits, etc},
90        u'Rigid bodies':{'data': {rigid body dict}},
91        u'Covariance':{'data':{covariance data dict}},
92        u'Controls':{'data':{controls data dict}},
93        u'Notebook':{'data':[notebook list]},
94        u'Restraints':{'data':{restraint data dict}},
95        u'Constraints':{'data':{constraint data dict}}]
96        }
97      nameList = [                #NB: reproduces tree order
98        [u'Notebook',],
99        [u'Controls',],
100        [u'Covariance',],
101        [u'Constraints',],
102        [u'Restraints',],
103        [u'Rigid bodies',],
104        [u'PWDR FAP.XRA Bank 1',
105             u'Comments',
106             u'Limits',
107             u'Background',
108             u'Instrument Parameters',
109             u'Sample Parameters',
110             u'Peak List',
111             u'Index Peak List',
112             u'Unit Cells List',
113             u'Reflection Lists'],
114        [u'Phases', u'fap']
115        ]
116    '''
117    # Let IOError be thrown if file does not exist
118    # if not ospath.exists(ProjFile):
119    #     print ('\n*** Error attempt to open project file that does not exist:\n   '+
120    #         str(ProjFile))
121    #     return
122    file = open(ProjFile,'rb')
123    # print('loading from file: {}'.format(ProjFile))
124    Project = {}
125    nameList = []
126    try:
127        while True:
128            try:
129                data = cPickle.load(file)
130            except EOFError:
131                break
132            datum = data[0]
133            Project[datum[0]] = {'data':datum[1]}
134            nameList.append([datum[0],])
135            for datus in data[1:]:
136                Project[datum[0]][datus[0]] = datus[1]
137                nameList[-1].append(datus[0])
138        file.close()
139        # print('project load successful')
140    except:
141        raise IOError("Error reading file "+str(ProjFile)+". This is not a GSAS-II .gpx file")
142    finally:
143        file.close()
144    return Project,nameList
145
146def SaveDictToProjFile(Project,nameList,ProjFile):
147    '''Save a GSAS-II project file from dictionary/nameList created by LoadDictFromProjFile
148
149    :param dict Project: representation of gpx file following the GSAS-II
150        tree structure as described for LoadDictFromProjFile
151    :param list nameList: names of main tree entries & subentries used to reconstruct project file
152    :param str ProjFile: full file name for output project.gpx file (including extension)
153    '''
154    file = open(ProjFile,'wb')
155    try:
156        for name in nameList:
157            data = []
158            item = Project[name[0]]
159            data.append([name[0],item['data']])
160            for item2 in name[1:]:
161                data.append([item2,item[item2]])
162            cPickle.dump(data,file,1)
163    finally:
164        file.close()
165
166def ImportPowder(reader,filename):
167    '''Use a reader to import a powder diffraction data file
168
169    :param str reader: a scriptable reader
170    :param str filename: full name of powder data file; can be "multi-Bank" data
171
172    :returns list rdlist: list of reader objects containing powder data, one for each
173        "Bank" of data encountered in file. Items in reader object of interest are:
174
175          * rd.comments: list of str: comments found on powder file
176          * rd.dnames: list of str: data nammes suitable for use in GSASII data tree NB: duplicated in all rd entries in rdlist
177          * rd.powderdata: list of numpy arrays: pos,int,wt,zeros,zeros,zeros as needed for a PWDR entry in  GSASII data tree.
178    '''
179    rdfile,rdpath,descr = imp.find_module(reader)
180    rdclass = imp.load_module(reader,rdfile,rdpath,descr)
181    rd = rdclass.GSAS_ReaderClass()
182    if not rd.scriptable:
183        print(u'**** ERROR: '+reader+u' is not a scriptable reader')
184        return None
185    fl = open(filename,'rb')
186    rdlist = []
187    if rd.ContentsValidator(fl):
188        fl.seek(0)
189        repeat = True
190        rdbuffer = {} # create temporary storage for file reader
191        block = 0
192        while repeat: # loop if the reader asks for another pass on the file
193            block += 1
194            repeat = False
195            rd.objname = ospath.basename(filename)
196            flag = rd.Reader(filename,fl,None,buffer=rdbuffer,blocknum=block,)
197            if flag:
198                rdlist.append(copy.deepcopy(rd)) # save the result before it is written over
199                if rd.repeat:
200                    repeat = True
201        return rdlist
202    print(rd.errors)
203    return None
204
205def SetDefaultDData(dType,histoName,NShkl=0,NDij=0):
206    '''Create an initial Histogram dictionary
207
208    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
209    '''
210    if dType in ['SXC','SNC']:
211        return {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
212            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
213            'Extinction':['Lorentzian','None', {'Tbar':0.1,'Cos2TM':0.955,
214            'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}],
215            'Flack':[0.0,False]}
216    elif dType == 'SNT':
217        return {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
218            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
219            'Extinction':['Lorentzian','None', {
220            'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}]}
221    elif 'P' in dType:
222        return {'Histogram':histoName,'Show':False,'Scale':[1.0,False],
223            'Pref.Ori.':['MD',1.0,False,[0,0,1],0,{},[],0.1],
224            'Size':['isotropic',[1.,1.,1.],[False,False,False],[0,0,1],
225                [1.,1.,1.,0.,0.,0.],6*[False,]],
226            'Mustrain':['isotropic',[1000.0,1000.0,1.0],[False,False,False],[0,0,1],
227                NShkl*[0.01,],NShkl*[False,]],
228            'HStrain':[NDij*[0.0,],NDij*[False,]],
229            'Extinction':[0.0,False],'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}}
230
231
232def PreSetup(data):
233    '''Create part of an initial (empty) phase dictionary
234
235    from GSASIIphsGUI.py, near end of UpdatePhaseData
236
237    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
238    '''
239    if 'RBModels' not in data:
240        data['RBModels'] = {}
241    if 'MCSA' not in data:
242        data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[.8,1.2],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}}
243    if 'dict' in str(type(data['MCSA']['Results'])):
244        data['MCSA']['Results'] = []
245    if 'Modulated' not in data['General']:
246        data['General']['Modulated'] = False
247    if 'modulated' in data['General']['Type']:
248        data['General']['Modulated'] = True
249        data['General']['Type'] = 'nuclear'
250
251
252def SetupGeneral(data, dirname):
253    """Helps initialize phase data.
254
255    From GSASIIphsGui.py, function of the same name. Minor changes for imports etc.
256
257    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
258    """
259    mapDefault = {'MapType':'','RefList':'','Resolution':0.5,'Show bonds':True,
260                'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.,'Flip':False}
261    generalData = data['General']
262    atomData = data['Atoms']
263    generalData['AtomTypes'] = []
264    generalData['Isotopes'] = {}
265
266    if 'Isotope' not in generalData:
267        generalData['Isotope'] = {}
268    if 'Data plot type' not in generalData:
269        generalData['Data plot type'] = 'Mustrain'
270    if 'POhkl' not in generalData:
271        generalData['POhkl'] = [0,0,1]
272    if 'Map' not in generalData:
273        generalData['Map'] = mapDefault.copy()
274    if 'Flip' not in generalData:
275        generalData['Flip'] = {'RefList':'','Resolution':0.5,'Norm element':'None',
276            'k-factor':0.1,'k-Max':20.,}
277    if 'testHKL' not in generalData['Flip']:
278        generalData['Flip']['testHKL'] = [[0,0,2],[2,0,0],[1,1,1],[0,2,0],[1,2,3]]
279    if 'doPawley' not in generalData:
280        generalData['doPawley'] = False     #ToDo: change to ''
281    if 'Pawley dmin' not in generalData:
282        generalData['Pawley dmin'] = 1.0
283    if 'Pawley neg wt' not in generalData:
284        generalData['Pawley neg wt'] = 0.0
285    if 'Algolrithm' in generalData.get('MCSA controls',{}) or \
286        'MCSA controls' not in generalData:
287        generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50],
288        'dmin':2.0,'Algorithm':'log','Jump coeff':[0.95,0.5],'boltzmann':1.0,
289        'fast parms':[1.0,1.0,1.0],'log slope':0.9,'Cycles':1,'Results':[],'newDmin':True}
290    if 'AtomPtrs' not in generalData:
291        generalData['AtomPtrs'] = [3,1,7,9]
292        if generalData['Type'] == 'macromolecular':
293            generalData['AtomPtrs'] = [6,4,10,12]
294        elif generalData['Type'] == 'magnetic':
295            generalData['AtomPtrs'] = [3,1,10,12]
296    if generalData['Type'] in ['modulated',]:
297        generalData['Modulated'] = True
298        generalData['Type'] = 'nuclear'
299        if 'Super' not in generalData:
300            generalData['Super'] = 1
301            generalData['SuperVec'] = [[0,0,.1],False,4]
302            generalData['SSGData'] = {}
303        if '4DmapData' not in generalData:
304            generalData['4DmapData'] = mapDefault.copy()
305            generalData['4DmapData'].update({'MapType':'Fobs'})
306    if 'Modulated' not in generalData:
307        generalData['Modulated'] = False
308    if 'HydIds' not in generalData:
309        generalData['HydIds'] = {}
310    cx,ct,cs,cia = generalData['AtomPtrs']
311    generalData['NoAtoms'] = {}
312    generalData['BondRadii'] = []
313    generalData['AngleRadii'] = []
314    generalData['vdWRadii'] = []
315    generalData['AtomMass'] = []
316    generalData['Color'] = []
317    if generalData['Type'] == 'magnetic':
318        generalData['MagDmin'] = generalData.get('MagDmin',1.0)
319        landeg = generalData.get('Lande g',[])
320    generalData['Mydir'] = dirname
321    badList = {}
322    for iat,atom in enumerate(atomData):
323        atom[ct] = atom[ct].lower().capitalize()              #force to standard form
324        if generalData['AtomTypes'].count(atom[ct]):
325            generalData['NoAtoms'][atom[ct]] += atom[cx+3]*float(atom[cs+1])
326        elif atom[ct] != 'UNK':
327            Info = G2elem.GetAtomInfo(atom[ct])
328            if not Info:
329                if atom[ct] not in badList:
330                    badList[atom[ct]] = 0
331                badList[atom[ct]] += 1
332                atom[ct] = 'UNK'
333                continue
334            atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo
335            generalData['AtomTypes'].append(atom[ct])
336            generalData['Z'] = Info['Z']
337            generalData['Isotopes'][atom[ct]] = Info['Isotopes']
338            generalData['BondRadii'].append(Info['Drad'])
339            generalData['AngleRadii'].append(Info['Arad'])
340            generalData['vdWRadii'].append(Info['Vdrad'])
341            if atom[ct] in generalData['Isotope']:
342                if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]:
343                    isotope = generalData['Isotopes'][atom[ct]].keys()[-1]
344                    generalData['Isotope'][atom[ct]] = isotope
345                generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][atom[ct]]]['Mass'])
346            else:
347                generalData['Isotope'][atom[ct]] = 'Nat. Abund.'
348                if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]:
349                    isotope = generalData['Isotopes'][atom[ct]].keys()[-1]
350                    generalData['Isotope'][atom[ct]] = isotope
351                generalData['AtomMass'].append(Info['Mass'])
352            generalData['NoAtoms'][atom[ct]] = atom[cx+3]*float(atom[cs+1])
353            generalData['Color'].append(Info['Color'])
354            if generalData['Type'] == 'magnetic':
355                if len(landeg) < len(generalData['AtomTypes']):
356                    landeg.append(2.0)
357    if generalData['Type'] == 'magnetic':
358        generalData['Lande g'] = landeg[:len(generalData['AtomTypes'])]
359
360    if badList:
361        msg = 'Warning: element symbol(s) not found:'
362        for key in badList:
363            msg += '\n\t' + key
364            if badList[key] > 1:
365                msg += ' (' + str(badList[key]) + ' times)'
366        raise Exception("Phase error:\n" + msg)
367        # wx.MessageBox(msg,caption='Element symbol error')
368    F000X = 0.
369    F000N = 0.
370    for i,elem in enumerate(generalData['AtomTypes']):
371        F000X += generalData['NoAtoms'][elem]*generalData['Z']
372        isotope = generalData['Isotope'][elem]
373        F000N += generalData['NoAtoms'][elem]*generalData['Isotopes'][elem][isotope]['SL'][0]
374    generalData['F000X'] = F000X
375    generalData['F000N'] = F000N
376    import GSASIImath as G2mth
377    generalData['Mass'] = G2mth.getMass(generalData)
378
379
380############################################
381############ CIF export helpers ############
382############################################
383## These functions are translated from
384## exports/G2export_CIF.py
385## They are used in G2Phase.export_CIF
386
387
388
389def make_empty_project(author=None, filename=None):
390    """Creates an dictionary in the style of GSASIIscriptable, for an empty
391    project.
392
393    If no author name or filename is supplied, 'no name' and
394    <current dir>/test_output.gpx are used , respectively.
395
396    Returns: project dictionary, name list
397
398    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
399    """
400    if not filename:
401        filename = os.path.join(os.getcwd(), 'test_output.gpx')
402    else:
403        filename = os.path.abspath(filename)
404    gsasii_version = str(GSASIIpath.GetVersionNumber())
405    LoadG2fil()
406    import matplotlib as mpl
407    python_library_versions = G2fil.get_python_versions([mpl, np, sp])
408
409    controls_data = dict(G2obj.DefaultControls)
410    controls_data['LastSavedAs'] = unicode(filename)
411    controls_data['LastSavedUsing'] = gsasii_version
412    controls_data['PythonVersions'] = python_library_versions
413    if author:
414        controls_data['Author'] = author
415
416    output = {'Constraints': {'data': {'HAP': [], 'Hist': [], 'Phase': [],
417                                       'Global': []}},
418              'Controls': {'data': controls_data},
419              u'Covariance': {'data': {}},
420              u'Notebook': {'data': ['']},
421              u'Restraints': {'data': {}},
422              u'Rigid bodies': {'data': {'RBIds': {'Residue': [], 'Vector': []},
423                                'Residue': {'AtInfo': {}},
424                                'Vector':  {'AtInfo': {}}}}}
425
426    names = [[u'Notebook'], [u'Controls'], [u'Covariance'],
427             [u'Constraints'], [u'Restraints'], [u'Rigid bodies']]
428
429    return output, names
430
431
432class G2ImportException(Exception):
433    pass
434
435
436def import_generic(filename, readerlist):
437    """Attempt to import a filename, using a list of reader objects.
438
439    Returns the first reader object which worked."""
440    # Translated from OnImportGeneric method in GSASII.py
441    primaryReaders, secondaryReaders = [], []
442    for reader in readerlist:
443        flag = reader.ExtensionValidator(filename)
444        if flag is None:
445            secondaryReaders.append(reader)
446        elif flag:
447            primaryReaders.append(reader)
448    if not secondaryReaders and not primaryReaders:
449        raise G2ImportException("Could not read file: ", filename)
450
451    with open(filename, 'Ur') as fp:
452        rd_list = []
453
454        for rd in primaryReaders + secondaryReaders:
455            # Initialize reader
456            rd.selections = []
457            rd.dnames = []
458            rd.ReInitialize()
459            # Rewind file
460            fp.seek(0)
461            rd.errors = ""
462            if not rd.ContentsValidator(fp):
463                # Report error
464                pass
465            if len(rd.selections) > 1:
466                # Select data?
467                # GSASII.py:543
468                raise G2ImportException("Not sure what data to select")
469
470            block = 0
471            rdbuffer = {}
472            fp.seek(0)
473            repeat = True
474            while repeat:
475                repeat = False
476                block += 1
477                rd.objname = os.path.basename(filename)
478                flag = rd.Reader(filename, fp, buffer=rdbuffer, blocknum=block)
479                if flag:
480                    # Omitting image loading special cases
481                    rd.readfilename = filename
482                    rd_list.append(copy.deepcopy(rd))
483                    repeat = rd.repeat
484                else:
485                    raise G2ImportException("{}.Reader() returned:".format(rd),
486                                            flag)
487
488            if rd_list:
489                if rd.warnings:
490                    print("Read warning by", rd.formatName, "reader:",
491                          rd.warnings, file=sys.stderr)
492                return rd_list
493    raise G2ImportException("No reader could read file: " + filename)
494
495
496def load_iprms(instfile, reader):
497    """Loads instrument parameters from a file, and edits the
498    given reader.
499
500    Returns a 2-tuple of (Iparm1, Iparm2) parameters
501    """
502    LoadG2fil()
503    ext = os.path.splitext(instfile)[1]
504
505    if ext.lower() == '.instprm':
506        # New GSAS File, load appropriately
507        with open(instfile) as f:
508            lines = f.readlines()
509        bank = reader.powderentry[2]
510        numbanks = reader.numbanks
511        iparms = G2fil.ReadPowderInstprm(lines, bank, numbanks, reader)
512
513        reader.instfile = instfile
514        reader.instmsg = 'GSAS-II file' + instfile
515        return iparms
516    elif ext.lower() not in ('.prm', '.inst', '.ins'):
517        raise ValueError('Expected .prm file, found: ', instfile)
518
519    # It's an old GSAS file, load appropriately
520    Iparm = {}
521    with open(instfile, 'Ur') as fp:
522        for line in fp:
523            if '#' in line:
524                continue
525            Iparm[line[:12]] = line[12:-1]
526    ibanks = int(Iparm.get('INS   BANK  ', '1').strip())
527    if ibanks == 1:
528        reader.instbank = 1
529        reader.powderentry[2] = 1
530        reader.instfile = instfile
531        reader.instmsg = instfile + ' bank ' + str(reader.instbank)
532        return G2fil.SetPowderInstParms(Iparm, reader)
533    # TODO handle >1 banks
534    raise NotImplementedError("Check GSASIIfiles.py: ReadPowderInstprm")
535
536def load_pwd_from_reader(reader, instprm, existingnames=[]):
537    """Loads powder data from a reader object, and assembles it into a GSASII data tree.
538
539    :returns: (name, tree) - 2-tuple of the histogram name (str), and data
540
541    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
542    """
543    HistName = 'PWDR ' + G2obj.StripUnicode(reader.idstring, '_')
544    HistName = unicode(G2obj.MakeUniqueLabel(HistName, existingnames))
545
546    try:
547        Iparm1, Iparm2 = instprm
548    except ValueError:
549        Iparm1, Iparm2 = load_iprms(instprm, reader)
550
551    Ymin = np.min(reader.powderdata[1])
552    Ymax = np.max(reader.powderdata[1])
553    valuesdict = {'wtFactor': 1.0,
554                  'Dummy': False,
555                  'ranId': ran.randint(0, sys.maxint),
556                  'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax,
557                  'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax,
558                  'qPlot': False, 'dPlot': False, 'sqrtPlot': False,
559                  'Yminmax': [Ymin, Ymax]}
560    reader.Sample['ranId'] = valuesdict['ranId']
561
562    # Ending keys:
563    # [u'Reflection Lists',
564    #  u'Limits',
565    #  'data',
566    #  u'Index Peak List',
567    #  u'Comments',
568    #  u'Unit Cells List',
569    #  u'Sample Parameters',
570    #  u'Peak List',
571    #  u'Background',
572    #  u'Instrument Parameters']
573    Tmin = np.min(reader.powderdata[0])
574    Tmax = np.max(reader.powderdata[0])
575
576    default_background = [['chebyschev', False, 3, 1.0, 0.0, 0.0],
577                          {'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0, 'peaksList': []}]
578
579    output_dict = {u'Reflection Lists': {},
580                   u'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin, Tmax]]),
581                   u'data': [valuesdict, reader.powderdata, HistName],
582                   u'Index Peak List': [[], []],
583                   u'Comments': reader.comments,
584                   u'Unit Cells List': [],
585                   u'Sample Parameters': reader.Sample,
586                   u'Peak List': {'peaks': [], 'sigDict': {}},
587                   u'Background': reader.pwdparms.get('Background', default_background),
588                   u'Instrument Parameters': [Iparm1, Iparm2],
589                   }
590
591    names = [u'Comments',
592             u'Limits',
593             u'Background',
594             u'Instrument Parameters',
595             u'Sample Parameters',
596             u'Peak List',
597             u'Index Peak List',
598             u'Unit Cells List',
599             u'Reflection Lists']
600
601    # TODO controls?? GSASII.py:1664-7
602
603    return HistName, [HistName] + names, output_dict
604
605
606def _deep_copy_into(from_, into):
607    """Helper function for reloading .gpx file. See G2Project.reload()
608
609    :author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
610    """
611    if isinstance(from_, dict) and isinstance(into, dict):
612        combined_keys = set(from_.keys()).union(into.keys())
613        for key in combined_keys:
614            if key in from_ and key in into:
615                both_dicts = (isinstance(from_[key], dict)
616                              and isinstance(into[key], dict))
617                both_lists = (isinstance(from_[key], list)
618                              and isinstance(into[key], list))
619                if both_dicts or both_lists:
620                    _deep_copy_into(from_[key], into[key])
621                else:
622                    into[key] = from_[key]
623            elif key in from_:
624                into[key] = from_[key]
625            else:  # key in into
626                del into[key]
627    elif isinstance(from_, list) and isinstance(into, list):
628        if len(from_) == len(into):
629            for i in xrange(len(from_)):
630                both_dicts = (isinstance(from_[i], dict)
631                              and isinstance(into[i], dict))
632                both_lists = (isinstance(from_[i], list)
633                              and isinstance(into[i], list))
634                if both_dicts or both_lists:
635                    _deep_copy_into(from_[i], into[i])
636                else:
637                    into[i] = from_[i]
638        else:
639            into[:] = from_
640
641
642class G2ObjectWrapper(object):
643    """Base class for all GSAS-II object wrappers.
644
645    The underlying GSAS-II format can be accessed as `wrapper.data`. A number
646    of overrides are implemented so that the wrapper behaves like a dictionary.
647
648    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
649    """
650    def __init__(self, datadict):
651        self.data = datadict
652
653    def __getitem__(self, key):
654        return self.data[key]
655
656    def __setitem__(self, key, value):
657        self.data[key] = value
658
659    def __contains__(self, key):
660        return key in self.data
661
662    def get(self, k, d=None):
663        return self.data.get(k, d)
664
665    def keys(self):
666        return self.data.keys()
667
668    def values(self):
669        return self.data.values()
670
671    def items(self):
672        return self.data.items()
673
674
675class G2Project(G2ObjectWrapper):
676    """Represents an entire GSAS-II project."""
677    def __init__(self, gpxfile=None, author=None, filename=None):
678        """Loads a GSAS-II project from a specified filename.
679
680        :param str gpxfile: Existing .gpx file to be loaded. If nonexistent,
681            creates an empty project.
682        :param str author: Author's name. Optional.
683        :param str filename: The filename the project should be saved to in
684            the future. If both filename and gpxfile are present, the project is
685            loaded from the gpxfile, then set to  save to filename in the future"""
686        if gpxfile is None:
687            filename = os.path.abspath(os.path.expanduser(filename))
688            self.filename = filename
689            self.data, self.names = make_empty_project(author=author, filename=filename)
690        elif isinstance(gpxfile, str):
691            # TODO set author, filename
692            self.filename = os.path.abspath(os.path.expanduser(gpxfile))
693            self.data, self.names = LoadDictFromProjFile(gpxfile)
694            self.update_ids()
695            self.index_ids()
696        else:
697            raise ValueError("Not sure what to do with gpxfile")
698
699    @classmethod
700    def from_dict_and_names(cls, gpxdict, names, filename=None):
701        """Creates a :class:`G2Project` directly from
702        a dictionary and a list of names. If in doubt, do not use this.
703
704        :returns: a :class:`G2Project`
705        """
706        out = cls()
707        if filename:
708            filename = os.path.abspath(os.path.expanduser(filename))
709            out.filename = filename
710            gpxdict['Controls']['data']['LastSavedAs'] = filename
711        else:
712            try:
713                out.filename = gpxdict['Controls']['data']['LastSavedAs']
714            except KeyError:
715                out.filename = None
716        out.data = gpxdict
717        out.names = names
718
719    def save(self, filename=None):
720        """Saves the project, either to the current filename, or to a new file.
721
722        Updates self.filename if a new filename provided"""
723        # TODO update LastSavedUsing ?
724        if filename:
725            filename = os.path.abspath(os.path.expanduser(filename))
726            self.data['Controls']['data']['LastSavedAs'] = filename
727            self.filename = filename
728        elif not self.filename:
729            raise AttributeError("No file name to save to")
730        SaveDictToProjFile(self.data, self.names, self.filename)
731
732    def add_powder_histogram(self, datafile, iparams):
733        """Loads a powder data histogram into the project.
734
735        Automatically checks for an instrument parameter file, or one can be
736        provided.
737
738        :param str datafile: The powder data file to read, a filename.
739        :param str iparams: The instrument parameters file, a filename.
740
741        :returns: A :class:`G2PwdrData` object representing
742            the histogram
743        """
744        LoadG2fil()
745        datafile = os.path.abspath(os.path.expanduser(datafile))
746        iparams = os.path.abspath(os.path.expanduser(iparams))
747        pwdrreaders = import_generic(datafile, PwdrDataReaders)
748        histname, new_names, pwdrdata = load_pwd_from_reader(
749                                          pwdrreaders[0], iparams,
750                                          [h.name for h in self.histograms()])
751        if histname in self.data:
752            print("Warning - redefining histogram", histname)
753        else:
754            if self.names[-1][0] == 'Phases':
755                self.names.insert(-1, new_names)
756            else:
757                self.names.append(new_names)
758        self.data[histname] = pwdrdata
759        self.update_ids()
760        return self.histogram(histname)
761
762    def add_phase(self, phasefile, phasename=None, histograms=[]):
763        """Loads a phase into the project from a .cif file
764
765        :param str phasefile: The CIF file from which to import the phase.
766        :param str phasename: The name of the new phase, or None for the default
767        :param list histograms: The names of the histograms to associate with
768            this phase
769
770        :returns: A :class:`G2Phase` object representing the
771            new phase.
772        """
773        LoadG2fil()
774        histograms = [self.histogram(h).name for h in histograms]
775        phasefile = os.path.abspath(os.path.expanduser(phasefile))
776
777        # TODO handle multiple phases in a file
778        phasereaders = import_generic(phasefile, PhaseReaders)
779        phasereader = phasereaders[0]
780
781        phasename = phasename or phasereader.Phase['General']['Name']
782        phaseNameList = [p.name for p in self.phases()]
783        phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList)
784        phasereader.Phase['General']['Name'] = phasename
785
786        if 'Phases' not in self.data:
787            self.data[u'Phases'] = { 'data': None }
788        assert phasename not in self.data['Phases'], "phase names should be unique"
789        self.data['Phases'][phasename] = phasereader.Phase
790
791        if phasereader.Constraints:
792            Constraints = self.data['Constraints']
793            for i in phasereader.Constraints:
794                if isinstance(i, dict):
795                    if '_Explain' not in Constraints:
796                        Constraints['_Explain'] = {}
797                    Constraints['_Explain'].update(i)
798                else:
799                    Constraints['Phase'].append(i)
800
801        data = self.data['Phases'][phasename]
802        generalData = data['General']
803        SGData = generalData['SGData']
804        NShkl = len(G2spc.MustrainNames(SGData))
805        NDij = len(G2spc.HStrainNames(SGData))
806        Super = generalData.get('Super', 0)
807        if Super:
808            SuperVec = np.array(generalData['SuperVec'][0])
809        else:
810            SuperVec = []
811        UseList = data['Histograms']
812
813        for hist in histograms:
814            self.link_histogram_phase(hist, phasename)
815            # if histname.startswith('HKLF '):
816            #     raise NotImplementedError("Does not support HKLF yet")
817            # elif histname.startswith('PWDR '):
818            #     hist = self.histogram(histname)
819            #     hist['Reflection Lists'][generalData['Name']] = {}
820            #     UseList[histname] = SetDefaultDData('PWDR', histname, NShkl=NShkl, NDij=NDij)
821            #     for key, val in [('Use', True), ('LeBail', False),
822            #                      ('newLeBail', True),
823            #                      ('Babinet', {'BabA': [0.0, False],
824            #                                   'BabU': [0.0, False]})]:
825            #         if key not in UseList[histname]:
826            #             UseList[histname][key] = val
827            # else:
828            #     raise NotImplementedError("Unexpected histogram" + histname)
829
830        for obj in self.names:
831            if obj[0] == 'Phases':
832                phasenames = obj
833                break
834        else:
835            phasenames = [u'Phases']
836            self.names.append(phasenames)
837        phasenames.append(unicode(phasename))
838
839        # TODO should it be self.filename, not phasefile?
840        SetupGeneral(data, os.path.dirname(phasefile))
841        self.index_ids()
842
843        self.update_ids()
844        return self.phase(phasename)
845
846    def link_histogram_phase(self, histogram, phase):
847        """Associates a given histogram and phase.
848
849        .. seealso::
850
851            :func:`~G2Project.histogram`
852            :func:`~G2Project.phase`"""
853        hist = self.histogram(histogram)
854        phase = self.phase(phase)
855
856        generalData = phase['General']
857
858        if hist.name.startswith('HKLF '):
859            raise NotImplementedError("HKLF not yet supported")
860        elif hist.name.startswith('PWDR '):
861            hist['Reflection Lists'][generalData['Name']] = {}
862            UseList = phase['Histograms']
863            SGData = generalData['SGData']
864            NShkl = len(G2spc.MustrainNames(SGData))
865            NDij = len(G2spc.HStrainNames(SGData))
866            UseList[hist.name] = SetDefaultDData('PWDR', hist.name, NShkl=NShkl, NDij=NDij)
867            UseList[hist.name]['hId'] = hist.id
868            for key, val in [('Use', True), ('LeBail', False),
869                             ('newLeBail', True),
870                             ('Babinet', {'BabA': [0.0, False],
871                                          'BabU': [0.0, False]})]:
872                if key not in UseList[hist.name]:
873                    UseList[hist.name][key] = val
874        else:
875            raise RuntimeError("Unexpected histogram" + hist.name)
876
877
878    def reload(self):
879        """Reload self from self.filename"""
880        data, names = LoadDictFromProjFile(self.filename)
881        self.names = names
882        # Need to deep copy the new data file data into the current tree,
883        # so that any existing G2Phase, or G2PwdrData objects will still be
884        # valid
885        _deep_copy_into(from_=data, into=self.data)
886
887    def refine(self, newfile=None, printFile=None, makeBack=False):
888        # TODO migrate to RefineCore
889        # G2strMain.RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
890        #      calcControls,pawleyLookup,ifPrint,printFile,dlg)
891        # index_ids will automatically save the project
892        self.index_ids()
893        # TODO G2strMain does not properly use printFile
894        G2strMain.Refine(self.filename, makeBack=makeBack)
895        # Reload yourself
896        self.reload()
897
898    def histogram(self, histname):
899        """Returns the histogram named histname, or None if it does not exist.
900
901        :param histname: The name of the histogram (str), or ranId or index.
902        :returns: A :class:`G2PwdrData` object, or None if
903            the histogram does not exist
904
905        .. seealso::
906            :func:`~GSASIIscriptable.G2Project.histograms`
907            :func:`~GSASIIscriptable.G2Project.phase`
908            :func:`~GSASIIscriptable.G2Project.phases`
909            """
910        if isinstance(histname, G2PwdrData):
911            if histname.proj == self:
912                return histname
913        if histname in self.data:
914            return G2PwdrData(self.data[histname], self)
915        for key, val in G2obj.HistIdLookup.items():
916            name, ranId = val
917            # histname can be either ranId (key) or index (val)
918            if ranId == histname or key == str(histname):
919                return self.histogram(name)
920
921    def histograms(self):
922        """Return a list of all histograms, as
923        :class:`G2PwdrData` objects
924
925        .. seealso::
926            :func:`~GSASIIscriptable.G2Project.histograms`
927            :func:`~GSASIIscriptable.G2Project.phase`
928            :func:`~GSASIIscriptable.G2Project.phases`
929            """
930        output = []
931        for obj in self.names:
932            if len(obj) > 1 and obj[0] != u'Phases':
933                output.append(self.histogram(obj[0]))
934        return output
935
936    def phase(self, phasename):
937        """
938        Gives an object representing the specified phase in this project.
939
940        :param str phasename: The name of the desired phase. Either the name
941            (str), the phase's ranId, or the phase's index
942        :returns: A :class:`G2Phase` object
943        :raises: KeyError
944
945        .. seealso::
946            :func:`~GSASIIscriptable.G2Project.histograms`
947            :func:`~GSASIIscriptable.G2Project.phase`
948            :func:`~GSASIIscriptable.G2Project.phases`
949            """
950        phases = self.data['Phases']
951        if phasename in phases:
952            return G2Phase(phases[phasename], phasename, self)
953        for key, val in G2obj.PhaseIdLookup.items():
954            name, ranId = val
955            # phasename can be either ranId (key) or index (val)
956            if ranId == phasename or key == str(phasename):
957                return self.phase(name)
958
959    def phases(self):
960        """
961        Returns a list of all the phases in the project.
962
963        :returns: A list of :class:`G2Phase` objects
964
965        .. seealso::
966            :func:`~GSASIIscriptable.G2Project.histogram`
967            :func:`~GSASIIscriptable.G2Project.histograms`
968            :func:`~GSASIIscriptable.G2Project.phase`
969            """
970        for obj in self.names:
971            if obj[0] == 'Phases':
972                return [self.phase(p) for p in obj[1:]]
973        return []
974
975    def update_ids(self):
976        """Makes sure all phases and histograms have proper hId and pId"""
977        # Translated from GetUsedHistogramsAndPhasesfromTree,
978        #   GSASIIdataGUI.py:4107
979        for i, h in enumerate(self.histograms()):
980            h.id = i
981        for i, p in enumerate(self.phases()):
982            p.id = i
983
984    def do_refinements(self, refinements, histogram='all', phase='all',
985                       outputnames=None, makeBack=False):
986        """Conducts a series of refinements. Wrapper around iter_refinements
987
988        :param list refinements: A list of dictionaries defining refinements
989        :param str histogram: Name of histogram for refinements to be applied
990            to, or 'all'
991        :param str phase: Name of phase for refinements to be applied to, or
992            'all'
993        """
994        for proj in self.iter_refinements(refinements, histogram, phase,
995                                          outputnames, makeBack):
996            pass
997        return self
998
999    def iter_refinements(self, refinements, histogram='all', phase='all',
1000                         outputnames=None, makeBack=False):
1001        """Conducts a series of refinements, iteratively. Stops after every
1002        refinement and yields this project, to allow error checking or
1003        logging of intermediate results.
1004
1005        >>> def checked_refinements(proj):
1006        ...     for p in proj.iter_refinements(refs):
1007        ...         # Track intermediate results
1008        ...         log(p.histogram('0').residuals)
1009        ...         log(p.phase('0').get_cell())
1010        ...         # Check if parameter diverged, nonsense answer, or whatever
1011        ...         if is_something_wrong(p):
1012        ...             raise Exception("I need a human!")
1013
1014        :param list refinements: A list of dictionaries defining refinements
1015        :param str histogram: Name of histogram for refinements to be applied
1016            to, or 'all'
1017        :param str phase: Name of phase for refinements to be applied to, or
1018            'all'
1019        """
1020        if outputnames:
1021            if len(refinements) != len(outputnames):
1022                raise ValueError("Should have same number of outuputs to"
1023                                 "refinements")
1024        else:
1025            outputnames = [None for r in refinements]
1026
1027        for output, refinement in zip(outputnames, refinements):
1028            self.set_refinement(refinement, histogram)
1029            # Handle 'once' args - refinements that are disabled after this
1030            # refinement
1031            if 'once' in refinement:
1032                temp = {'set': refinement['once']}
1033                self.set_refinement(temp, histogram, phase)
1034
1035            if output:
1036                self.save(output)
1037
1038            self.refine(makeBack=makeBack)
1039            yield self
1040
1041            # Handle 'once' args - refinements that are disabled after this
1042            # refinement
1043            if 'once' in refinement:
1044                temp = {'clear': refinement['once']}
1045                self.set_refinement(temp, histogram, phase)
1046
1047    def set_refinement(self, refinement, histogram='all', phase='all'):
1048        """Apply specified refinements to a given histogram(s) or phase(s).
1049
1050        Refinement parameters are categorize in three groups:
1051
1052        1. Histogram parameters
1053        2. Phase parameters
1054        3. Histogram-and-Phase (HAP) parameters
1055
1056        :param dict refinement: The refinements to be conducted
1057        :param histogram: Either a name of a histogram (str), a list of
1058            histogram names, or 'all' (default)
1059        :param phase: Either a name of a phase (str), a list of phase names, or
1060            'all' (default)
1061
1062        .. seealso::
1063            :func:`~G2PwdrData.set_refinements`
1064            :func:`~G2PwdrData.clear_refinements`
1065            :func:`~G2Phase.set_refinements`
1066            :func:`~G2Phase.clear_refinements`
1067            :func:`~G2Phase.set_HAP_refinements`
1068            :func:`~G2Phase.clear_HAP_refinements`"""
1069
1070        if histogram == 'all':
1071            hists = self.histograms()
1072        elif isinstance(histogram, str) or isinstance(histogram, int):
1073            hists = [self.histogram(histogram)]
1074        else:
1075            hists = [self.histogram(name) for name in histogram]
1076
1077        if phase == 'all':
1078            phases = self.phases()
1079        elif isinstance(phase, str) or isinstance(phase, int):
1080            phases = [self.phase(phase)]
1081        else:
1082            phases = [self.phase(name) for name in phase]
1083
1084
1085        # TODO: HAP parameters:
1086        #   Babinet
1087        #   Extinction
1088        #   HStrain
1089        #   Mustrain
1090        #   Pref. Ori
1091        #   Size
1092
1093        pwdr_set = {}
1094        phase_set = {}
1095        hap_set = {}
1096        for key, val in refinement.get('set', {}).items():
1097            # Apply refinement options
1098            if G2PwdrData.is_valid_refinement_key(key):
1099                pwdr_set[key] = val
1100            elif G2Phase.is_valid_refinement_key(key):
1101                phase_set[key] = val
1102            elif G2Phase.is_valid_HAP_refinement_key(key):
1103                hap_set[key] = val
1104            else:
1105                raise ValueError("Unknown refinement key", key)
1106
1107        for hist in hists:
1108            hist.set_refinements(pwdr_set)
1109        for phase in phases:
1110            phase.set_refinements(phase_set)
1111        for phase in phases:
1112            phase.set_HAP_refinements(hap_set, hists)
1113
1114        pwdr_clear = {}
1115        phase_clear = {}
1116        hap_clear = {}
1117        for key, val in refinement.get('clear', {}).items():
1118            # Clear refinement options
1119            if G2PwdrData.is_valid_refinement_key(key):
1120                pwdr_clear[key] = val
1121            elif G2Phase.is_valid_refinement_key(key):
1122                phase_clear[key] = val
1123            elif G2Phase.is_valid_HAP_refinement_key(key):
1124                hap_set[key] = val
1125            else:
1126                raise ValueError("Unknown refinement key", key)
1127
1128        for hist in hists:
1129            hist.clear_refinements(pwdr_clear)
1130        for phase in phases:
1131            phase.clear_refinements(phase_clear)
1132        for phase in phases:
1133            phase.clear_HAP_refinements(hap_clear, hists)
1134
1135    def index_ids(self):
1136        import GSASIIstrIO as G2strIO
1137        self.save()
1138        return G2strIO.GetUsedHistogramsAndPhases(self.filename)
1139
1140    def add_constraint_raw(self, cons_scope, constr):
1141        """Adds a constraint of type consType to the project.
1142        cons_scope should be one of "Hist", "Phase", "HAP", or "Global".
1143
1144        WARNING it does not check the constraint is well-constructed"""
1145        constrs = self.data['Constraints']['data']
1146        if 'Global' not in constrs:
1147            constrs['Global'] = []
1148        constrs[cons_scope].append(constr)
1149
1150    def hold_many(self, vars, type):
1151        """Apply holds for all the variables in vars, for constraint of a given type.
1152
1153        type is passed directly to add_constraint_raw as consType
1154
1155        :param list vars: A list of variables to hold. Either :class:`GSASIIobj.G2VarObj` objects,
1156            string variable specifiers, or arguments for :meth:`make_var_obj`
1157        :param str type: A string constraint type specifier. See
1158        :class:`~GSASIIscriptable.G2Project.add_constraint_raw`"""
1159        for var in vars:
1160            if isinstance(var, str):
1161                var = self.make_var_obj(var)
1162            elif not isinstance(var, G2obj.G2VarObj):
1163                var = self.make_var_obj(*var)
1164            self.add_constraint_raw(type, [[1.0, var], None, None, 'h'])
1165
1166    def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None,
1167                     reloadIdx=True):
1168        """Wrapper to create a G2VarObj. Takes either a string representaiton ("p:h:name:a")
1169        or individual names of phase, histogram, varname, and atomId.
1170
1171        Automatically converts string phase, hist, or atom names into the ID required
1172        by G2VarObj."""
1173
1174        if reloadIdx:
1175            self.index_ids()
1176
1177        # If string representation, short circuit
1178        if hist is None and varname is None and atomId is None:
1179            if isinstance(phase, str) and ':' in phase:
1180                return G2obj.G2VarObj(phase)
1181
1182        # Get phase index
1183        phaseObj = None
1184        if isinstance(phase, G2Phase):
1185            phaseObj = phase
1186            phase = G2obj.PhaseRanIdLookup[phase.ranId]
1187        elif phase in self.data['Phases']:
1188            phaseObj = self.phase(phase)
1189            phaseRanId = phaseObj.ranId
1190            phase = G2obj.PhaseRanIdLookup[phaseRanId]
1191
1192        # Get histogram index
1193        if isinstance(hist, G2PwdrData):
1194            hist = G2obj.HistRanIdLookup[hist.ranId]
1195        elif hist in self.data:
1196            histRanId = self.histogram(hist).ranId
1197            hist = G2obj.HistRanIdLookup[histRanId]
1198
1199        # Get atom index (if any)
1200        if isinstance(atomId, G2AtomRecord):
1201            atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId]
1202        elif phaseObj:
1203            atomObj = phaseObj.atom(atomId)
1204            if atomObj:
1205                atomRanId = atomObj.ranId
1206                atomId = G2obj.AtomRanIdLookup[phase][atomRanId]
1207
1208        return G2obj.G2VarObj(phase, hist, varname, atomId)
1209
1210
1211class G2AtomRecord(G2ObjectWrapper):
1212    """Wrapper for an atom record. Has convenient accessors via @property.
1213
1214
1215    Available accessors: label, type, refinement_flags, coordinates,
1216        occupancy, ranId, id, adp_flag, uiso
1217
1218    Example:
1219
1220    >>> atom = some_phase.atom("O3")
1221    >>> # We can access the underlying data format
1222    >>> atom.data
1223    ['O3', 'O-2', '', ... ]
1224    >>> # We can also use wrapper accessors
1225    >>> atom.coordinates
1226    (0.33, 0.15, 0.5)
1227    >>> atom.refinement_flags
1228    u'FX'
1229    >>> atom.ranId
1230    4615973324315876477
1231    >>> atom.occupancy
1232    1.0
1233    """
1234    def __init__(self, data, indices, proj):
1235        self.data = data
1236        self.cx, self.ct, self.cs, self.cia = indices
1237        self.proj = proj
1238
1239    @property
1240    def label(self):
1241        return self.data[self.ct-1]
1242
1243    @property
1244    def type(self):
1245        return self.data[self.ct]
1246
1247    @property
1248    def refinement_flags(self):
1249        return self.data[self.ct+1]
1250
1251    @refinement_flags.setter
1252    def refinement_flags(self, other):
1253        # Automatically check it is a valid refinement
1254        for c in other:
1255            if c not in ' FXU':
1256                raise ValueError("Invalid atom refinement: ", other)
1257        self.data[self.ct+1] = unicode(other)
1258
1259    @property
1260    def coordinates(self):
1261        return tuple(self.data[self.cx:self.cx+3])
1262
1263    @property
1264    def occupancy(self):
1265        return self.data[self.cx+3]
1266
1267    @occupancy.setter
1268    def occupancy(self, val):
1269        self.data[self.cx+3] = float(val)
1270
1271    @property
1272    def ranId(self):
1273        return self.data[self.cia+8]
1274
1275    @property
1276    def adp_flag(self):
1277        # Either 'I' or 'A'
1278        return self.data[self.cia]
1279
1280    @property
1281    def uiso(self):
1282        if self.adp_flag == 'I':
1283            return self.data[self.cia+1]
1284        else:
1285            return self.data[self.cia+2:self.cia+8]
1286
1287
1288class G2PwdrData(G2ObjectWrapper):
1289    """Wraps a Powder Data Histogram."""
1290    def __init__(self, data, proj):
1291        self.data = data
1292        self.proj = proj
1293
1294    @staticmethod
1295    def is_valid_refinement_key(key):
1296        valid_keys = ['Limits', 'Sample Parameters', 'Background',
1297                      'Instrument Parameters']
1298        return key in valid_keys
1299
1300    @property
1301    def name(self):
1302        return self.data['data'][-1]
1303
1304    @property
1305    def ranId(self):
1306        return self.data['data'][0]['ranId']
1307
1308    @property
1309    def residuals(self):
1310        data = self.data['data'][0]
1311        return {key: data[key]
1312                for key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']}
1313
1314    @property
1315    def id(self):
1316        self.proj.update_ids()
1317        return self.data['data'][0]['hId']
1318
1319    @id.setter
1320    def id(self, val):
1321        self.data['data'][0]['hId'] = val
1322
1323    def fit_fixed_points(self):
1324        """Attempts to apply a background fit to the fixed points currently specified."""
1325        def SetInstParms(Inst):
1326            dataType = Inst['Type'][0]
1327            insVary = []
1328            insNames = []
1329            insVals = []
1330            for parm in Inst:
1331                insNames.append(parm)
1332                insVals.append(Inst[parm][1])
1333                if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1334                    'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1335                        Inst[parm][2] = False
1336            instDict = dict(zip(insNames, insVals))
1337            instDict['X'] = max(instDict['X'], 0.01)
1338            instDict['Y'] = max(instDict['Y'], 0.01)
1339            if 'SH/L' in instDict:
1340                instDict['SH/L'] = max(instDict['SH/L'], 0.002)
1341            return dataType, instDict, insVary
1342
1343        bgrnd = self.data['Background']
1344
1345        # Need our fixed points in order
1346        bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
1347        X = [x for x, y in bgrnd[1]['FixedPoints']]
1348        Y = [y for x, y in bgrnd[1]['FixedPoints']]
1349
1350        limits = self.data['Limits'][1]
1351        if X[0] > limits[0]:
1352            X = [limits[0]] + X
1353            Y = [Y[0]] + Y
1354        if X[-1] < limits[1]:
1355            X += [limits[1]]
1356            Y += [Y[-1]]
1357
1358        # Some simple lookups
1359        controls = self.proj['Controls']['data']
1360        inst, inst2 = self.data['Instrument Parameters']
1361        pwddata = self.data['data'][1]
1362
1363        # Construct the data for background fitting
1364        xBeg = np.searchsorted(pwddata[0], limits[0])
1365        xFin = np.searchsorted(pwddata[0], limits[1])
1366        xdata = pwddata[0][xBeg:xFin]
1367        ydata = si.interp1d(X,Y)(ma.getdata(xdata))
1368
1369        W = [1]*len(xdata)
1370        Z = [0]*len(xdata)
1371
1372        dataType, insDict, insVary = SetInstParms(inst)
1373        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1374
1375        # Do the fit
1376        data = np.array([xdata, ydata, W, Z, Z, Z])
1377        G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
1378                        prevVaryList=bakVary, controls=controls)
1379
1380        # Post-fit
1381        parmDict = {}
1382        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1383        parmDict.update(bakDict)
1384        parmDict.update(insDict)
1385        pwddata[3][xBeg:xFin] *= 0
1386        pwddata[5][xBeg:xFin] *= 0
1387        pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
1388
1389        # TODO adjust pwddata? GSASIIpwdGUI.py:1041
1390        # TODO update background
1391        self.proj.save()
1392
1393    def y_calc(self):
1394        return self.data['data'][1][3]
1395
1396    def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True):
1397        try:
1398            import matplotlib.pyplot as plt
1399            data = self.data['data'][1]
1400            if Yobs:
1401                plt.plot(data[0], data[1], label='Yobs')
1402            if Ycalc:
1403                plt.plot(data[0], data[3], label='Ycalc')
1404            if Background:
1405                plt.plot(data[0], data[4], label='Background')
1406            if Residual:
1407                plt.plot(data[0], data[5], label="Residual")
1408        except ImportError:
1409            pass
1410
1411    def set_refinements(self, refs):
1412        """Sets the refinement parameter 'key' to the specification 'value'
1413
1414        :param dict refs: A dictionary of the parameters to be set.
1415
1416        :returns: None"""
1417        do_fit_fixed_points = False
1418        for key, value in refs.items():
1419            if key == 'Limits':
1420                old_limits = self.data['Limits'][1]
1421                new_limits = value
1422                if isinstance(new_limits, dict):
1423                    if 'low' in new_limits:
1424                        old_limits[0] = new_limits['low']
1425                    if 'high' in new_limits:
1426                        old_limits[1] = new_limits['high']
1427                else:
1428                    old_limits[0], old_limits[1] = new_limits
1429            elif key == 'Sample Parameters':
1430                sample = self.data['Sample Parameters']
1431                for sparam in value:
1432                    sample[sparam][1] = True
1433            elif key == 'Background':
1434                bkg, peaks = self.data['Background']
1435
1436                # If True or False, just set the refine parameter
1437                if value in (True, False):
1438                    bkg[1] = value
1439                    return
1440
1441                if 'type' in value:
1442                    bkg[0] = value['type']
1443                if 'refine' in value:
1444                    bkg[1] = value['refine']
1445                if 'no. coeffs' in value:
1446                    cur_coeffs = bkg[2]
1447                    n_coeffs = value['no. coeffs']
1448                    if n_coeffs > cur_coeffs:
1449                        for x in range(n_coeffs - cur_coeffs):
1450                            bkg.append(0.0)
1451                    else:
1452                        for _ in range(cur_coeffs - n_coeffs):
1453                            bkg.pop()
1454                    bkg[2] = n_coeffs
1455                if 'coeffs' in value:
1456                    bkg[3:] = value['coeffs']
1457                if 'FixedPoints' in value:
1458                    peaks['FixedPoints'] = [(float(a), float(b))
1459                                            for a, b in value['FixedPoints']]
1460                if value.get('fit fixed points', False):
1461                    do_fit_fixed_points = True
1462
1463            elif key == 'Instrument Parameters':
1464                instrument, secondary = self.data['Instrument Parameters']
1465                for iparam in value:
1466                    try:
1467                        instrument[iparam][2] = True
1468                    except IndexError:
1469                        raise ValueError("Invalid key:", iparam)
1470            else:
1471                raise ValueError("Unknown key:", key)
1472        # Fit fixed points after the fact - ensure they are after fixed points
1473        # are added
1474        if do_fit_fixed_points:
1475            # Background won't be fit if refinement flag not set
1476            orig = self.data['Background'][0][1]
1477            self.data['Background'][0][1] = True
1478            self.fit_fixed_points()
1479            # Restore the previous value
1480            self.data['Background'][0][1] = orig
1481
1482    def clear_refinements(self, refs):
1483        """Clears the refinement parameter 'key' and its associated value.
1484
1485        :param dict refs: A dictionary of parameters to clear."""
1486        for key, value in refs.items():
1487            if key == 'Limits':
1488                old_limits, cur_limits = self.data['Limits']
1489                cur_limits[0], cur_limits[1] = old_limits
1490            elif key == 'Sample Parameters':
1491                sample = self.data['Sample Parameters']
1492                for sparam in value:
1493                    sample[sparam][1] = False
1494            elif key == 'Background':
1495                bkg, peaks = self.data['Background']
1496
1497                # If True or False, just set the refine parameter
1498                if value in (True, False):
1499                    bkg[1] = False
1500                    return
1501
1502                bkg[1] = False
1503                if 'FixedPoints' in value:
1504                    if 'FixedPoints' in peaks:
1505                        del peaks['FixedPoints']
1506            elif key == 'Instrument Parameters':
1507                instrument, secondary = self.data['Instrument Parameters']
1508                for iparam in value:
1509                    instrument[iparam][2] = False
1510            else:
1511                raise ValueError("Unknown key:", key)
1512
1513
1514class G2Phase(G2ObjectWrapper):
1515    """A wrapper object around a given phase.
1516
1517    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1518    """
1519    def __init__(self, data, name, proj):
1520        self.data = data
1521        self.name = name
1522        self.proj = proj
1523
1524    @staticmethod
1525    def is_valid_refinement_key(key):
1526        valid_keys = ["Cell", "Atoms", "LeBail"]
1527        return key in valid_keys
1528
1529    @staticmethod
1530    def is_valid_HAP_refinement_key(key):
1531        valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain",
1532                      "Pref.Ori.", "Show", "Size", "Use", "Scale"]
1533        return key in valid_keys
1534
1535    def atom(self, atomlabel):
1536        """Returns the atom specified by atomlabel, or None if it does not
1537        exist.
1538
1539        :param str atomlabel: The name of the atom (e.g. "O2")
1540        :returns: A :class:`G2AtomRecord` object
1541            representing the atom.
1542        """
1543        # Consult GSASIIobj.py for the meaning of this
1544        cx, ct, cs, cia = self.data['General']['AtomPtrs']
1545        ptrs = [cx, ct, cs, cia]
1546        atoms = self.data['Atoms']
1547        for atom in atoms:
1548            if atom[ct-1] == atomlabel:
1549                return G2AtomRecord(atom, ptrs, self.proj)
1550
1551    def atoms(self):
1552        """Returns a list of atoms present in the phase.
1553
1554        :returns: A list of :class:`G2AtomRecord` objects.
1555
1556        .. seealso::
1557            :func:`~GSASIIscriptable.G2Phase.atom`
1558            :class:`G2AtomRecord`
1559        """
1560        ptrs = self.data['General']['AtomPtrs']
1561        output = []
1562        atoms = self.data['Atoms']
1563        for atom in atoms:
1564            output.append(G2AtomRecord(atom, ptrs, self.proj))
1565        return output
1566
1567    def histograms(self):
1568        output = []
1569        for hname in self.data.get('Histograms', {}).keys():
1570            output.append(self.proj.histogram(hname))
1571        return output
1572
1573    @property
1574    def ranId(self):
1575        return self.data['ranId']
1576
1577    @property
1578    def id(self):
1579        return self.data['pId']
1580
1581    @id.setter
1582    def id(self, val):
1583        self.data['pId'] = val
1584
1585    def get_cell(self):
1586        """Returns a dictionary of the cell parameters, with keys:
1587            'length_a', 'length_b', 'length_c', 'angle_alpha', 'angle_beta', 'angle_gamma', 'volume'
1588
1589        :returns: a dict
1590
1591        .. seealso::
1592           :func:`~G2Phase.get_cell_and_esd`
1593        """
1594        cell = self.data['General']['Cell']
1595        return {'length_a': cell[1], 'length_b': cell[2], 'length_c': cell[3],
1596                'angle_alpha': cell[4], 'angle_beta': cell[5], 'angle_gamma': cell[6],
1597                'volume': cell[7]}
1598
1599    def get_cell_and_esd(self):
1600        """
1601        Returns a pair of dictionaries, the first representing the unit cell, the second
1602        representing the estimated standard deviations of the unit cell.
1603
1604        :returns: a tuple of two dictionaries
1605        .. seealso::
1606           :func:`~G2Phase.get_cell`
1607        """
1608        # translated from GSASIIstrIO.ExportBaseclass.GetCell
1609        import GSASIIstrIO as G2stIO
1610        import GSASIIlattice as G2lat
1611        import GSASIImapvars as G2mv
1612        try:
1613            pfx = str(self.id) + '::'
1614            sgdata = self['General']['SGData']
1615            covDict = self.proj['Covariance']['data']
1616
1617            parmDict = dict(zip(covDict.get('varyList',[]),
1618                                covDict.get('variables',[])))
1619            sigDict = dict(zip(covDict.get('varyList',[]),
1620                               covDict.get('sig',[])))
1621
1622            if covDict.get('covMatrix') is not None:
1623                sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
1624                                                  covDict['varyList'],
1625                                                  parmDict))
1626
1627            A, sigA = G2stIO.cellFill(pfx, sgdata, parmDict, sigDict)
1628            cellSig = G2stIO.getCellEsd(pfx, sgdata, A, self.proj['Covariance']['data'])
1629            cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),)
1630            cellDict, cellSigDict = {}, {}
1631            for i, key in enumerate(['length_a', 'length_b', 'length_c',
1632                                     'angle_alpha', 'angle_beta', 'angle_gamma',
1633                                     'volume']):
1634                cellDict[key] = cellList[i]
1635                cellSigDict[key] = cellSig[i]
1636            return cellDict, cellSigDict
1637        except KeyError:
1638            cell = self.get_cell()
1639            return cell, {key: 0.0 for key in cell}
1640
1641    def export_CIF(self, outputname, quickmode=True):
1642        """Write this phase to a .cif file named outputname
1643
1644        :param str outputname: The name of the .cif file to write to
1645        :param bool quickmode: Currently ignored. Carryover from exports.G2export_CIF"""
1646        # This code is all taken from exports/G2export_CIF.py
1647        # Functions copied have the same names
1648        import GSASIImath as G2mth
1649        import GSASIImapvars as G2mv
1650        from exports import G2export_CIF as cif
1651
1652        CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M")
1653        CIFname = os.path.splitext(self.proj.filename)[0]
1654        CIFname = os.path.split(CIFname)[1]
1655        CIFname = ''.join([c if ord(c) < 128 else ''
1656                           for c in CIFname.replace(' ', '_')])
1657        try:
1658            author = self.proj['Controls']['data'].get('Author','').strip()
1659        except KeyError:
1660            pass
1661        oneblock = True
1662
1663        covDict = self.proj['Covariance']['data']
1664        parmDict = dict(zip(covDict.get('varyList',[]),
1665                            covDict.get('variables',[])))
1666        sigDict = dict(zip(covDict.get('varyList',[]),
1667                           covDict.get('sig',[])))
1668
1669        if covDict.get('covMatrix') is not None:
1670            sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
1671                                              covDict['varyList'],
1672                                              parmDict))
1673
1674        with open(outputname, 'w') as fp:
1675            fp.write(' \n' + 70*'#' + '\n')
1676            cif.WriteCIFitem(fp, 'data_' + CIFname)
1677            # from exports.G2export_CIF.WritePhaseInfo
1678            cif.WriteCIFitem(fp, '\n# phase info for '+str(self.name) + ' follows')
1679            cif.WriteCIFitem(fp, '_pd_phase_name', self.name)
1680            # TODO get esds
1681            cellDict = self.get_cell()
1682            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
1683            names = ['length_a','length_b','length_c',
1684                     'angle_alpha','angle_beta ','angle_gamma',
1685                     'volume']
1686            for key, val in cellDict.items():
1687                cif.WriteCIFitem(fp, '_cell_' + key, G2mth.ValEsd(val))
1688
1689            cif.WriteCIFitem(fp, '_symmetry_cell_setting',
1690                         self.data['General']['SGData']['SGSys'])
1691
1692            spacegroup = self.data['General']['SGData']['SpGrp'].strip()
1693            # regularize capitalization and remove trailing H/R
1694            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
1695            cif.WriteCIFitem(fp, '_symmetry_space_group_name_H-M', spacegroup)
1696
1697            # generate symmetry operations including centering and center of symmetry
1698            SymOpList, offsetList, symOpList, G2oprList, G2opcodes = G2spc.AllOps(
1699                self.data['General']['SGData'])
1700            cif.WriteCIFitem(fp, 'loop_\n    _space_group_symop_id\n    _space_group_symop_operation_xyz')
1701            for i, op in enumerate(SymOpList,start=1):
1702                cif.WriteCIFitem(fp, '   {:3d}  {:}'.format(i,op.lower()))
1703
1704            # TODO skipped histograms, exports/G2export_CIF.py:880
1705
1706            # report atom params
1707            if self.data['General']['Type'] in ['nuclear','macromolecular']:        #this needs macromolecular variant, etc!
1708                cif.WriteAtomsNuclear(fp, self.data, self.name, parmDict, sigDict, [])
1709                # self._WriteAtomsNuclear(fp, parmDict, sigDict)
1710            else:
1711                raise Exception,"no export for "+str(self.data['General']['Type'])+" coordinates implemented"
1712            # report cell contents
1713            cif.WriteComposition(fp, self.data, self.name, parmDict)
1714            if not quickmode and self.data['General']['Type'] == 'nuclear':      # report distances and angles
1715                # WriteDistances(fp,self.name,SymOpList,offsetList,symOpList,G2oprList)
1716                raise NotImplementedError("only quickmode currently supported")
1717            if 'Map' in self.data['General'] and 'minmax' in self.data['General']['Map']:
1718                cif.WriteCIFitem(fp,'\n# Difference density results')
1719                MinMax = self.data['General']['Map']['minmax']
1720                cif.WriteCIFitem(fp,'_refine_diff_density_max',G2mth.ValEsd(MinMax[0],-0.009))
1721                cif.WriteCIFitem(fp,'_refine_diff_density_min',G2mth.ValEsd(MinMax[1],-0.009))
1722
1723
1724    def set_refinements(self, refs):
1725        """Sets the refinement parameter 'key' to the specification 'value'
1726
1727        :param dict refs: A dictionary of the parameters to be set.
1728
1729        :returns: None"""
1730        for key, value in refs.items():
1731            if key == "Cell":
1732                self.data['General']['Cell'][0] = True
1733            elif key == "Atoms":
1734                cx, ct, cs, cia = self.data['General']['AtomPtrs']
1735
1736                for atomlabel, atomrefinement in value.items():
1737                    if atomlabel == 'all':
1738                        for atom in self.atoms():
1739                            atom.refinement_flags = atomrefinement
1740                    else:
1741                        atom = self.atom(atomlabel)
1742                        atom.refinement_flags = atomrefinement
1743            elif key == "LeBail":
1744                hists = self.data['Histograms']
1745                for hname, hoptions in hists.items():
1746                    if 'LeBail' not in hoptions:
1747                        hoptions['newLeBail'] = True
1748                    hoptions['LeBail'] = bool(value)
1749            else:
1750                raise ValueError("Unknown key:", key)
1751
1752    def clear_refinements(self, refs):
1753        """Clears a given set of parameters.
1754
1755        :param dict refs: The parameters to clear"""
1756        for key, value in refs.items():
1757            if key == "Cell":
1758                self.data['General']['Cell'][0] = False
1759            elif key == "Atoms":
1760                cx, ct, cs, cia = self.data['General']['AtomPtrs']
1761
1762                for atomlabel in value:
1763                    atom = self.atom(atomlabel)
1764                    # Set refinement to none
1765                    atom.refinement_flags = ' '
1766            elif key == "LeBail":
1767                hists = self.data['Histograms']
1768                for hname, hoptions in hists.items():
1769                    if 'LeBail' not in hoptions:
1770                        hoptions['newLeBail'] = True
1771                    hoptions['LeBail'] = False
1772            else:
1773                raise ValueError("Unknown key:", key)
1774
1775    def set_HAP_refinements(self, refs, histograms='all'):
1776        """Sets the given HAP refinement parameters between this phase and
1777        the given histograms
1778
1779        :param dict refs: A dictionary of the parameters to be set.
1780        :param histograms: Either 'all' (default) or a list of the histograms
1781            whose HAP parameters will be set with this phase. Histogram and phase
1782            must already be associated
1783
1784        :returns: None
1785        """
1786        if histograms == 'all':
1787            histograms = self.data['Histograms'].values()
1788        else:
1789            histograms = [self.data['Histograms'][h.name] for h in histograms]
1790
1791        for key, val in refs.items():
1792            for h in histograms:
1793                if key == 'Babinet':
1794                    try:
1795                        sets = list(val)
1796                    except ValueError:
1797                        sets = ['BabA', 'BabU']
1798                    for param in sets:
1799                        if param not in ['BabA', 'BabU']:
1800                            raise ValueError("Not sure what to do with" + param)
1801                        for hist in histograms:
1802                            hist['Babinet'][param][1] = True
1803                elif key == 'Extinction':
1804                    for h in histograms:
1805                        h['Extinction'][1] = True
1806                elif key == 'HStrain':
1807                    for h in histograms:
1808                        hist['HStrain'][1] = [True for p in hist['Hstrain'][0]]
1809                elif key == 'Mustrain':
1810                    for h in histograms:
1811                        mustrain = h['Mustrain']
1812                        newType = None
1813                        if isinstance(val, str):
1814                            if val in ['isotropic', 'uniaxial', 'generalized']:
1815                                newType = val
1816                            else:
1817                                raise ValueError("Not a Mustrain type: " + val)
1818                        elif isinstance(val, dict):
1819                            newType = val.get('type', None)
1820                            direction = val.get('direction', None)
1821
1822                        if newType:
1823                            mustrain[0] = newType
1824                            if newType == 'isotropic':
1825                                mustrain[2][0] = True
1826                                mustrain[5] = [False for p in mustrain[4]]
1827                            elif newType == 'uniaxial':
1828                                pass
1829                            else:  # newtype == 'generalized'
1830                                mustrain[2] = [False for p in mustrain[1]]
1831                        if direction:
1832                            direction = [int(n) for n in direction]
1833                            if len(direction) != 3:
1834                                raise ValueError("Expected hkl, found", direction)
1835                            mustrain[3] = direction
1836                elif key == 'Pref.Ori.':
1837                    for h in histograms:
1838                        h['Pref.Ori.'][2] = True
1839                elif key == 'Show':
1840                    for h in histograms:
1841                        h['Show'] = True
1842                elif key == 'Size':
1843                    raise NotImplementedError()
1844                elif key == 'Use':
1845                    for h in histograms:
1846                        h['Use'] = True
1847                elif key == 'Scale':
1848                    for h in histograms:
1849                        h['Scale'][1] = False
1850
1851    def clear_HAP_refinements(self, refs, histograms='all'):
1852        """Clears the given HAP refinement parameters between this phase and
1853        the given histograms
1854
1855        :param dict refs: A dictionary of the parameters to be cleared.
1856        :param histograms: Either 'all' (default) or a list of the histograms
1857            whose HAP parameters will be cleared with this phase. Histogram and
1858            phase must already be associated
1859
1860        :returns: None
1861        """
1862        if histograms == 'all':
1863            histograms = self.data['Histograms'].values()
1864        else:
1865            histograms = [self.data['Histograms'][h.name] for h in histograms]
1866
1867        for key, val in refs.items():
1868            for h in histograms:
1869                if key == 'Babinet':
1870                    try:
1871                        sets = list(val)
1872                    except ValueError:
1873                        sets = ['BabA', 'BabU']
1874                    for param in sets:
1875                        if param not in ['BabA', 'BabU']:
1876                            raise ValueError("Not sure what to do with" + param)
1877                        for hist in histograms:
1878                            hist['Babinet'][param][1] = False
1879                elif key == 'Extinction':
1880                    for h in histograms:
1881                        h['Extinction'][1] = False
1882                elif key == 'HStrain':
1883                    for h in histograms:
1884                        hist['HStrain'][1] = [False for p in hist['Hstrain'][0]]
1885                elif key == 'Mustrain':
1886                    for h in histograms:
1887                        mustrain = h['Mustrain']
1888                        mustrain[2] = [False for p in mustrain[1]]
1889                        mustrain[5] = [False for p in mustrain[4]]
1890                elif key == 'Pref.Ori.':
1891                    for h in histograms:
1892                        h['Pref.Ori.'][2] = False
1893                elif key == 'Show':
1894                    for h in histograms:
1895                        h['Show'] = False
1896                elif key == 'Size':
1897                    raise NotImplementedError()
1898                elif key == 'Use':
1899                    for h in histograms:
1900                        h['Use'] = False
1901                elif key == 'Scale':
1902                    for h in histograms:
1903                        h['Scale'][1] = False
1904
1905
1906##########################
1907# Command Line Interface #
1908##########################
1909
1910
1911def create(args):
1912    """The create subcommand."""
1913    proj = G2Project(filename=args.filename)
1914
1915    hist_objs = []
1916    for h in args.histograms:
1917        hist_objs.append(proj.add_powder_histogram(h, args.iparams))
1918
1919    for p in args.phases:
1920        proj.add_phase(p, histograms=hist_objs)
1921
1922    proj.save()
1923
1924
1925def dump(args):
1926    """The dump subcommand"""
1927    if not args.histograms and not args.phases:
1928        args.raw = True
1929    if args.raw:
1930        import IPython.lib.pretty as pretty
1931
1932    for fname in args.files:
1933        if args.raw:
1934            proj, nameList = LoadDictFromProjFile(fname)
1935            print("file:", fname)
1936            print("names:", nameList)
1937            for key, val in proj.items():
1938                print(key, ":")
1939                pretty.pprint(val)
1940        else:
1941            proj = G2Project(fname)
1942            if args.histograms:
1943                hists = proj.histograms()
1944                for h in hists:
1945                    print(fname, "hist", h.id, h.name)
1946            if args.phases:
1947                phase_list = proj.phases()
1948                for p in phase_list:
1949                    print(fname, "phase", p.id, p.name)
1950
1951
1952def IPyBrowse(args):
1953    """Load a .gpx file and then open a IPython shell to browse it
1954    """
1955    for fname in args.files:
1956        proj, nameList = LoadDictFromProjFile(fname)
1957        msg = "\nfname {} loaded into proj (dict) with names in nameList".format(fname)
1958        GSASIIpath.IPyBreak_base(msg)
1959        break
1960
1961
1962def refine(args):
1963    """The refine subcommand"""
1964    proj = G2Project(args.gpxfile)
1965    if args.refinements is None:
1966        proj.refine()
1967    else:
1968        import json
1969        with open(args.refinements) as refs:
1970            refs = json.load(refs)
1971        proj.do_refinements(refs['refinements'])
1972
1973
1974def seqrefine(args):
1975    """The seqrefine subcommand"""
1976    raise NotImplementedError("seqrefine is not yet implemented")
1977
1978
1979def export(args):
1980    """The export subcommand"""
1981    # Export phase as CIF to args.exportfile
1982    proj = G2Project(args.gpxfile)
1983    phase = proj.phase(args.phase)
1984    phase.export_CIF(args.exportfile)
1985
1986
1987def _args_kwargs(*args, **kwargs):
1988    return args, kwargs
1989
1990# A dictionary of the name of each subcommand, and a tuple
1991# of its associated function and a list of its arguments
1992# The arguments are passed directly to the add_argument() method
1993# of an argparse.ArgumentParser
1994subcommands = {"create":
1995               (create, [_args_kwargs('filename',
1996                                      help='the project file to create. should end in .gpx'),
1997
1998                         _args_kwargs('-i', '--iparams',
1999                                      help='instrument parameter file'),
2000
2001                         _args_kwargs('-d', '--histograms',
2002                                      nargs='+',
2003                                      help='list of datafiles to add as histograms'),
2004
2005                         _args_kwargs('-p', '--phases',
2006                                      nargs='+',
2007                                      help='list of phases to add. phases are '
2008                                           'automatically associated with all '
2009                                           'histograms given.')]),
2010
2011               "dump": (dump, [_args_kwargs('-d', '--histograms',
2012                                     action='store_true',
2013                                     help='list histograms in files, overrides --raw'),
2014
2015                               _args_kwargs('-p', '--phases',
2016                                            action='store_true',
2017                                            help='list phases in files, overrides --raw'),
2018
2019                               _args_kwargs('-r', '--raw',
2020                                      action='store_true', help='dump raw file contents, default'),
2021
2022                               _args_kwargs('files', nargs='+')]),
2023
2024               "refine":
2025               (refine, [_args_kwargs('gpxfile', help='the project file to refine'),
2026                         _args_kwargs('refinements',
2027                                      help='json file of refinements to apply. if not present'
2028                                           ' refines file as-is',
2029                                      default=None,
2030                                      nargs='?')]),
2031
2032               "seqrefine": (seqrefine, []),
2033               "export": (export, [_args_kwargs('gpxfile',
2034                                                help='the project file from which to export'),
2035                                   _args_kwargs('phase', help='identifier of phase to export'),
2036                                   _args_kwargs('exportfile', help='the .cif file to export to')]),
2037               "browse": (IPyBrowse, [_args_kwargs('files', nargs='+',
2038                                                   help='list of files to browse')])}
2039
2040
2041def main():
2042    '''The command line interface for GSASIIscriptable.
2043
2044    Executes one of the following subcommands:
2045
2046        * :func:`create`
2047        * :func:`dump`
2048        * :func:`refine`
2049        * :func:`seqrefine`
2050        * :func:`export`
2051        * :func:`browse`
2052
2053    .. seealso::
2054        :func:`create`
2055        :func:`dump`
2056        :func:`refine`
2057        :func:`seqrefine`
2058        :func:`export`
2059        :func:`browse`
2060    '''
2061    parser = argparse.ArgumentParser()
2062    subs = parser.add_subparsers()
2063
2064    # Create all of the specified subparsers
2065    for name, (func, args) in subcommands.items():
2066        new_parser = subs.add_parser(name)
2067        for listargs, kwds in args:
2068            new_parser.add_argument(*listargs, **kwds)
2069        new_parser.set_defaults(func=func)
2070
2071    # Parse and trigger subcommand
2072    result = parser.parse_args()
2073    result.func(result)
2074
2075    # argv = sys.argv
2076    # if len(argv) > 1 and argv[1] in subcommands:
2077    #     subcommands[argv[1]](*argv[2:])
2078    # elif len(argv) == 1 or argv[1] in ('help', '--help', '-h'):
2079    #     # TODO print usage
2080    #     subcommand_names = ' | '.join(sorted(subcommands.keys()))
2081    #     print("USAGE: {} [ {} ] ...".format(argv[0], subcommand_names))
2082    # else:
2083    #     print("Unknown subcommand: {}".format(argv[1]))
2084    #     print("Available subcommands:")
2085    #     for name in sorted(subcommands.keys()):
2086    #         print("\t{}".format(name))
2087    #     sys.exit(-1)
2088    # sys.exit(0)
2089
2090if __name__ == '__main__':
2091    main()
2092
2093
2094    # from refinements.py
2095#      USAGE = """USAGE: {} datafile instparams phasefile projectname refinements
2096
2097# datafile:    Input powder data
2098# intparams:   Corresponding instrument parameter file
2099# phasefile:   Phase to refine against data
2100# projectname: Project file to be created, should end in .gpx
2101# refinements: JSON file of refinements to be executed
2102# """
2103#     try:
2104#         datafile, instprm, phasefile, projectname, refinements = sys.argv[1:]
2105#     except ValueError:
2106#         print(USAGE.format(sys.argv[0]))
2107#         sys.exit(-1)
2108
2109#     try:
2110#         with open(refinements) as f:
2111#             refinements = json.load(f)
2112#     except IOError:
2113#         print("No such refinements file: {}".format(refinements))
2114
2115#     print("Creating project file \"{}\"...".format(projectname))
2116#     proj = G2Project(filename=projectname)
2117#     # Add the histogram
2118#     hist = proj.add_powder_histogram(datafile, instprm)
2119#     # Add the phase, and associate it with the histogram
2120#     proj.add_phase(phasefile, histograms=[hist.name])
2121
2122#     proj.do_refinements(refinements['refinements'])
2123#     proj.save()
2124
2125
2126    # from gpx_dumper
2127    # import IPython.lib.pretty as pretty
2128    # proj, nameList = LoadDictFromProjFile(sys.argv[1])
2129    # print("names:", nameList)
2130    # for key, val in proj.items():
2131    #     print(key, ":", sep='')
2132    #     pretty.pprint(val)
Note: See TracBrowser for help on using the repository browser.