source: trunk/GSASIIscriptable.py @ 3091

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

more documentation for G2Project

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