source: trunk/GSASIIscriptable.py @ 3085

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

removed comments, started documenting refinement specifiers

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