source: trunk/GSASIIscriptable.py @ 3098

Last change on this file since 3098 was 3098, checked in by toby, 5 years ago

fix powder data import errors with G2scriptable
File size: 91.3 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3########### SVN repository information ###################
4# $Date: 2017-04-12 15:12:45 -0500 (Wed, 12 Apr 2017) $
5# $Author: vondreele $
6# $Revision: 2777 $
7# $URL: https://subversion.xray.aps.anl.gov/pyGSAS/trunk/GSASIIscriptable.py $
8# $Id: GSASIIIO.py 2777 2017-04-12 20:12:45Z vondreele $
9########### SVN repository information ###################
10"""
11*GSASIIscriptable: Scripting Tools*
12-----------------------------------
13
14Routines for reading, writing, modifying and creating GSAS-II project (.gpx) files.
15
16Supports a command line interface as well. 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, True to refine
187Show                                         Boolean, True to refine
188Size                                         Not implemented
189Use                                          Boolean, True to refine
190Scale                                        Phase fraction; Boolean, True 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                try:
634                    flag = rd.Reader(filename, fp, buffer=rdbuffer, blocknum=block)
635                except:
636                    flag = False
637                if flag:
638                    # Omitting image loading special cases
639                    rd.readfilename = filename
640                    rd_list.append(copy.deepcopy(rd))
641                    repeat = rd.repeat
642                else:
643                    if GSASIIpath.GetConfigValue('debug'): print("{} Reader failed to read {}".format(rd.formatName,filename))
644            if rd_list:
645                if rd.warnings:
646                    print("Read warning by", rd.formatName, "reader:",
647                          rd.warnings, file=sys.stderr)
648                else:
649                    print("{} read by Reader {}\n".format(filename,rd.formatName))                   
650                return rd_list
651    raise G2ImportException("No reader could read file: " + filename)
652
653
654def load_iprms(instfile, reader):
655    """Loads instrument parameters from a file, and edits the
656    given reader.
657
658    Returns a 2-tuple of (Iparm1, Iparm2) parameters
659    """
660    LoadG2fil()
661    ext = os.path.splitext(instfile)[1]
662
663    if ext.lower() == '.instprm':
664        # New GSAS File, load appropriately
665        with open(instfile) as f:
666            lines = f.readlines()
667        bank = reader.powderentry[2]
668        numbanks = reader.numbanks
669        iparms = G2fil.ReadPowderInstprm(lines, bank, numbanks, reader)
670
671        reader.instfile = instfile
672        reader.instmsg = 'GSAS-II file' + instfile
673        return iparms
674    elif ext.lower() not in ('.prm', '.inst', '.ins'):
675        raise ValueError('Expected .prm file, found: ', instfile)
676
677    # It's an old GSAS file, load appropriately
678    Iparm = {}
679    with open(instfile, 'Ur') as fp:
680        for line in fp:
681            if '#' in line:
682                continue
683            Iparm[line[:12]] = line[12:-1]
684    ibanks = int(Iparm.get('INS   BANK  ', '1').strip())
685    if ibanks == 1:
686        reader.instbank = 1
687        reader.powderentry[2] = 1
688        reader.instfile = instfile
689        reader.instmsg = instfile + ' bank ' + str(reader.instbank)
690        return G2fil.SetPowderInstParms(Iparm, reader)
691    # TODO handle >1 banks
692    raise NotImplementedError("Check GSASIIfiles.py: ReadPowderInstprm")
693
694def load_pwd_from_reader(reader, instprm, existingnames=[]):
695    """Loads powder data from a reader object, and assembles it into a GSASII data tree.
696
697    :returns: (name, tree) - 2-tuple of the histogram name (str), and data
698
699    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
700    """
701    HistName = 'PWDR ' + G2obj.StripUnicode(reader.idstring, '_')
702    HistName = unicode(G2obj.MakeUniqueLabel(HistName, existingnames))
703
704    try:
705        Iparm1, Iparm2 = instprm
706    except ValueError:
707        Iparm1, Iparm2 = load_iprms(instprm, reader)
708
709    Ymin = np.min(reader.powderdata[1])
710    Ymax = np.max(reader.powderdata[1])
711    valuesdict = {'wtFactor': 1.0,
712                  'Dummy': False,
713                  'ranId': ran.randint(0, sys.maxint),
714                  'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax,
715                  'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax,
716                  'qPlot': False, 'dPlot': False, 'sqrtPlot': False,
717                  'Yminmax': [Ymin, Ymax]}
718    reader.Sample['ranId'] = valuesdict['ranId']
719
720    # Ending keys:
721    # [u'Reflection Lists',
722    #  u'Limits',
723    #  'data',
724    #  u'Index Peak List',
725    #  u'Comments',
726    #  u'Unit Cells List',
727    #  u'Sample Parameters',
728    #  u'Peak List',
729    #  u'Background',
730    #  u'Instrument Parameters']
731    Tmin = np.min(reader.powderdata[0])
732    Tmax = np.max(reader.powderdata[0])
733
734    default_background = [['chebyschev', False, 3, 1.0, 0.0, 0.0],
735                          {'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0, 'peaksList': []}]
736
737    output_dict = {u'Reflection Lists': {},
738                   u'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin, Tmax]]),
739                   u'data': [valuesdict, reader.powderdata, HistName],
740                   u'Index Peak List': [[], []],
741                   u'Comments': reader.comments,
742                   u'Unit Cells List': [],
743                   u'Sample Parameters': reader.Sample,
744                   u'Peak List': {'peaks': [], 'sigDict': {}},
745                   u'Background': reader.pwdparms.get('Background', default_background),
746                   u'Instrument Parameters': [Iparm1, Iparm2],
747                   }
748
749    names = [u'Comments',
750             u'Limits',
751             u'Background',
752             u'Instrument Parameters',
753             u'Sample Parameters',
754             u'Peak List',
755             u'Index Peak List',
756             u'Unit Cells List',
757             u'Reflection Lists']
758
759    # TODO controls?? GSASII.py:1664-7
760
761    return HistName, [HistName] + names, output_dict
762
763
764def _deep_copy_into(from_, into):
765    """Helper function for reloading .gpx file. See G2Project.reload()
766
767    :author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
768    """
769    if isinstance(from_, dict) and isinstance(into, dict):
770        combined_keys = set(from_.keys()).union(into.keys())
771        for key in combined_keys:
772            if key in from_ and key in into:
773                both_dicts = (isinstance(from_[key], dict)
774                              and isinstance(into[key], dict))
775                both_lists = (isinstance(from_[key], list)
776                              and isinstance(into[key], list))
777                if both_dicts or both_lists:
778                    _deep_copy_into(from_[key], into[key])
779                else:
780                    into[key] = from_[key]
781            elif key in from_:
782                into[key] = from_[key]
783            else:  # key in into
784                del into[key]
785    elif isinstance(from_, list) and isinstance(into, list):
786        if len(from_) == len(into):
787            for i in xrange(len(from_)):
788                both_dicts = (isinstance(from_[i], dict)
789                              and isinstance(into[i], dict))
790                both_lists = (isinstance(from_[i], list)
791                              and isinstance(into[i], list))
792                if both_dicts or both_lists:
793                    _deep_copy_into(from_[i], into[i])
794                else:
795                    into[i] = from_[i]
796        else:
797            into[:] = from_
798
799
800class G2ObjectWrapper(object):
801    """Base class for all GSAS-II object wrappers.
802
803    The underlying GSAS-II format can be accessed as `wrapper.data`. A number
804    of overrides are implemented so that the wrapper behaves like a dictionary.
805
806    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
807    """
808    def __init__(self, datadict):
809        self.data = datadict
810
811    def __getitem__(self, key):
812        return self.data[key]
813
814    def __setitem__(self, key, value):
815        self.data[key] = value
816
817    def __contains__(self, key):
818        return key in self.data
819
820    def get(self, k, d=None):
821        return self.data.get(k, d)
822
823    def keys(self):
824        return self.data.keys()
825
826    def values(self):
827        return self.data.values()
828
829    def items(self):
830        return self.data.items()
831
832
833class G2Project(G2ObjectWrapper):
834    """
835    Represents an entire GSAS-II project.
836
837    There are two ways to initialize it:
838
839    >>> # Load an existing project file
840    >>> proj = G2Project('filename.gpx')
841    >>> # Create a new project
842    >>> proj = G2Project(filename='new_file.gpx')
843    >>> # Specify an author
844    >>> proj = G2Project(author='Dr. So-And-So', filename='my_data.gpx')
845
846    Histograms can be accessed easily.
847
848    >>> # By name
849    >>> hist = proj.histogram('PWDR my-histogram-name')
850    >>> # Or by index
851    >>> hist = proj.histogram(0)
852    >>> assert hist.id == 0
853    >>> # Or by random id
854    >>> assert hist == proj.histogram(hist.ranId)
855
856    Phases can be accessed the same way.
857
858    >>> phase = proj.phase('name of phase')
859
860    New data can also be loaded via :func:`~G2Project.add_phase` and
861    :func:`~G2Project.add_powder_histogram`.
862
863    >>> hist = proj.add_powder_histogram('some_data_file.chi',
864                                         'instrument_parameters.prm')
865    >>> phase = proj.add_phase('my_phase.cif', histograms=[hist])
866
867    Parameters for Rietveld refinement can be turned on and off as well.
868    See :func:`~G2Project.set_refinement`, :func:`~G2Project.refine`,
869    :func:`~G2Project.iter_refinements`, :func:`~G2Project.do_refinements`.
870    """
871    def __init__(self, gpxfile=None, author=None, filename=None):
872        """Loads a GSAS-II project from a specified filename.
873
874        :param str gpxfile: Existing .gpx file to be loaded. If nonexistent,
875            creates an empty project.
876        :param str author: Author's name. Optional.
877        :param str filename: The filename the project should be saved to in
878            the future. If both filename and gpxfile are present, the project is
879            loaded from the gpxfile, then set to  save to filename in the future"""
880        if gpxfile is None:
881            filename = os.path.abspath(os.path.expanduser(filename))
882            self.filename = filename
883            self.data, self.names = make_empty_project(author=author, filename=filename)
884        elif isinstance(gpxfile, str):
885            # TODO set author, filename
886            self.filename = os.path.abspath(os.path.expanduser(gpxfile))
887            self.data, self.names = LoadDictFromProjFile(gpxfile)
888            self.update_ids()
889        else:
890            raise ValueError("Not sure what to do with gpxfile")
891
892    @classmethod
893    def from_dict_and_names(cls, gpxdict, names, filename=None):
894        """Creates a :class:`G2Project` directly from
895        a dictionary and a list of names. If in doubt, do not use this.
896
897        :returns: a :class:`G2Project`
898        """
899        out = cls()
900        if filename:
901            filename = os.path.abspath(os.path.expanduser(filename))
902            out.filename = filename
903            gpxdict['Controls']['data']['LastSavedAs'] = filename
904        else:
905            try:
906                out.filename = gpxdict['Controls']['data']['LastSavedAs']
907            except KeyError:
908                out.filename = None
909        out.data = gpxdict
910        out.names = names
911
912    def save(self, filename=None):
913        """Saves the project, either to the current filename, or to a new file.
914
915        Updates self.filename if a new filename provided"""
916        # TODO update LastSavedUsing ?
917        if filename:
918            filename = os.path.abspath(os.path.expanduser(filename))
919            self.data['Controls']['data']['LastSavedAs'] = filename
920            self.filename = filename
921        elif not self.filename:
922            raise AttributeError("No file name to save to")
923        SaveDictToProjFile(self.data, self.names, self.filename)
924
925    def add_powder_histogram(self, datafile, iparams, phases=[]):
926        """Loads a powder data histogram into the project.
927
928        Automatically checks for an instrument parameter file, or one can be
929        provided. Note that in unix fashion, "~" can be used to indicate the
930        home directory (e.g. ~/G2data/data.fxye).
931
932        :param str datafile: The powder data file to read, a filename.
933        :param str iparams: The instrument parameters file, a filename.
934        :param list phases: Phases to link to the new histogram
935
936        :returns: A :class:`G2PwdrData` object representing
937            the histogram
938        """
939        LoadG2fil()
940        datafile = os.path.abspath(os.path.expanduser(datafile))
941        iparams = os.path.abspath(os.path.expanduser(iparams))
942        pwdrreaders = import_generic(datafile, PwdrDataReaders)
943        histname, new_names, pwdrdata = load_pwd_from_reader(
944                                          pwdrreaders[0], iparams,
945                                          [h.name for h in self.histograms()])
946        if histname in self.data:
947            print("Warning - redefining histogram", histname)
948        elif self.names[-1][0] == 'Phases':
949            self.names.insert(-1, new_names)
950        else:
951            self.names.append(new_names)
952        self.data[histname] = pwdrdata
953        self.update_ids()
954
955        for phase in phases:
956            phase = self.phase(phase)
957            self.link_histogram_phase(histname, phase)
958
959        return self.histogram(histname)
960
961    def add_phase(self, phasefile, phasename=None, histograms=[]):
962        """Loads a phase into the project from a .cif file
963
964        :param str phasefile: The CIF file from which to import the phase.
965        :param str phasename: The name of the new phase, or None for the default
966        :param list histograms: The names of the histograms to associate with
967            this phase
968
969        :returns: A :class:`G2Phase` object representing the
970            new phase.
971        """
972        LoadG2fil()
973        histograms = [self.histogram(h).name for h in histograms]
974        phasefile = os.path.abspath(os.path.expanduser(phasefile))
975
976        # TODO handle multiple phases in a file
977        phasereaders = import_generic(phasefile, PhaseReaders)
978        phasereader = phasereaders[0]
979       
980        phasename = phasename or phasereader.Phase['General']['Name']
981        phaseNameList = [p.name for p in self.phases()]
982        phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList)
983        phasereader.Phase['General']['Name'] = phasename
984
985        if 'Phases' not in self.data:
986            self.data[u'Phases'] = { 'data': None }
987        assert phasename not in self.data['Phases'], "phase names should be unique"
988        self.data['Phases'][phasename] = phasereader.Phase
989
990        if phasereader.Constraints:
991            Constraints = self.data['Constraints']
992            for i in phasereader.Constraints:
993                if isinstance(i, dict):
994                    if '_Explain' not in Constraints:
995                        Constraints['_Explain'] = {}
996                    Constraints['_Explain'].update(i)
997                else:
998                    Constraints['Phase'].append(i)
999
1000        data = self.data['Phases'][phasename]
1001        generalData = data['General']
1002        SGData = generalData['SGData']
1003        NShkl = len(G2spc.MustrainNames(SGData))
1004        NDij = len(G2spc.HStrainNames(SGData))
1005        Super = generalData.get('Super', 0)
1006        if Super:
1007            SuperVec = np.array(generalData['SuperVec'][0])
1008        else:
1009            SuperVec = []
1010        UseList = data['Histograms']
1011
1012        for hist in histograms:
1013            self.link_histogram_phase(hist, phasename)
1014
1015        for obj in self.names:
1016            if obj[0] == 'Phases':
1017                phasenames = obj
1018                break
1019        else:
1020            phasenames = [u'Phases']
1021            self.names.append(phasenames)
1022        phasenames.append(unicode(phasename))
1023
1024        # TODO should it be self.filename, not phasefile?
1025        SetupGeneral(data, os.path.dirname(phasefile))
1026        self.index_ids()
1027
1028        self.update_ids()
1029        return self.phase(phasename)
1030
1031    def link_histogram_phase(self, histogram, phase):
1032        """Associates a given histogram and phase.
1033
1034        .. seealso::
1035
1036            :func:`~G2Project.histogram`
1037            :func:`~G2Project.phase`"""
1038        hist = self.histogram(histogram)
1039        phase = self.phase(phase)
1040
1041        generalData = phase['General']
1042
1043        if hist.name.startswith('HKLF '):
1044            raise NotImplementedError("HKLF not yet supported")
1045        elif hist.name.startswith('PWDR '):
1046            hist['Reflection Lists'][generalData['Name']] = {}
1047            UseList = phase['Histograms']
1048            SGData = generalData['SGData']
1049            NShkl = len(G2spc.MustrainNames(SGData))
1050            NDij = len(G2spc.HStrainNames(SGData))
1051            UseList[hist.name] = SetDefaultDData('PWDR', hist.name, NShkl=NShkl, NDij=NDij)
1052            UseList[hist.name]['hId'] = hist.id
1053            for key, val in [('Use', True), ('LeBail', False),
1054                             ('newLeBail', True),
1055                             ('Babinet', {'BabA': [0.0, False],
1056                                          'BabU': [0.0, False]})]:
1057                if key not in UseList[hist.name]:
1058                    UseList[hist.name][key] = val
1059        else:
1060            raise RuntimeError("Unexpected histogram" + hist.name)
1061
1062
1063    def reload(self):
1064        """Reload self from self.filename"""
1065        data, names = LoadDictFromProjFile(self.filename)
1066        self.names = names
1067        # Need to deep copy the new data file data into the current tree,
1068        # so that any existing G2Phase, or G2PwdrData objects will still be
1069        # valid
1070        _deep_copy_into(from_=data, into=self.data)
1071
1072    def refine(self, newfile=None, printFile=None, makeBack=False):
1073        # TODO migrate to RefineCore
1074        # G2strMain.RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
1075        #      calcControls,pawleyLookup,ifPrint,printFile,dlg)
1076        # index_ids will automatically save the project
1077        self.index_ids()
1078        # TODO G2strMain does not properly use printFile
1079        G2strMain.Refine(self.filename, makeBack=makeBack)
1080        # Reload yourself
1081        self.reload()
1082
1083    def histogram(self, histname):
1084        """Returns the histogram named histname, or None if it does not exist.
1085
1086        :param histname: The name of the histogram (str), or ranId or index.
1087        :returns: A :class:`G2PwdrData` object, or None if
1088            the histogram does not exist
1089
1090        .. seealso::
1091            :func:`~GSASIIscriptable.G2Project.histograms`
1092            :func:`~GSASIIscriptable.G2Project.phase`
1093            :func:`~GSASIIscriptable.G2Project.phases`
1094            """
1095        if isinstance(histname, G2PwdrData):
1096            if histname.proj == self:
1097                return histname
1098        if histname in self.data:
1099            return G2PwdrData(self.data[histname], self)
1100        try:
1101            # see if histname is an id or ranId
1102            histname = int(histname)
1103        except ValueError:
1104            return
1105
1106        for histogram in self.histograms():
1107            if histogram.id == histname or histogram.ranId == histname:
1108                return histogram
1109
1110    def histograms(self):
1111        """Return a list of all histograms, as
1112        :class:`G2PwdrData` objects
1113
1114        .. seealso::
1115            :func:`~GSASIIscriptable.G2Project.histograms`
1116            :func:`~GSASIIscriptable.G2Project.phase`
1117            :func:`~GSASIIscriptable.G2Project.phases`
1118            """
1119        output = []
1120        for obj in self.names:
1121            if len(obj) > 1 and obj[0] != u'Phases':
1122                output.append(self.histogram(obj[0]))
1123        return output
1124
1125    def phase(self, phasename):
1126        """
1127        Gives an object representing the specified phase in this project.
1128
1129        :param str phasename: The name of the desired phase. Either the name
1130            (str), the phase's ranId, or the phase's index
1131        :returns: A :class:`G2Phase` object
1132        :raises: KeyError
1133
1134        .. seealso::
1135            :func:`~GSASIIscriptable.G2Project.histograms`
1136            :func:`~GSASIIscriptable.G2Project.phase`
1137            :func:`~GSASIIscriptable.G2Project.phases`
1138            """
1139        phases = self.data['Phases']
1140        if phasename in phases:
1141            return G2Phase(phases[phasename], phasename, self)
1142
1143        try:
1144            # phasename should be phase index or ranId
1145            phasename = int(phasename)
1146        except ValueError:
1147            return
1148
1149        for phase in self.phases():
1150            if phase.id == phasename or phase.ranId == phasename:
1151                return phase
1152
1153    def phases(self):
1154        """
1155        Returns a list of all the phases in the project.
1156
1157        :returns: A list of :class:`G2Phase` objects
1158
1159        .. seealso::
1160            :func:`~GSASIIscriptable.G2Project.histogram`
1161            :func:`~GSASIIscriptable.G2Project.histograms`
1162            :func:`~GSASIIscriptable.G2Project.phase`
1163            """
1164        for obj in self.names:
1165            if obj[0] == 'Phases':
1166                return [self.phase(p) for p in obj[1:]]
1167        return []
1168
1169    def update_ids(self):
1170        """Makes sure all phases and histograms have proper hId and pId"""
1171        # Translated from GetUsedHistogramsAndPhasesfromTree,
1172        #   GSASIIdataGUI.py:4107
1173        for i, h in enumerate(self.histograms()):
1174            h.id = i
1175        for i, p in enumerate(self.phases()):
1176            p.id = i
1177
1178    def do_refinements(self, refinements, histogram='all', phase='all',
1179                       outputnames=None, makeBack=False):
1180        """Conducts a series of refinements. Wrapper around iter_refinements
1181
1182        :param list refinements: A list of dictionaries defining refinements
1183        :param str histogram: Name of histogram for refinements to be applied
1184            to, or 'all'
1185        :param str phase: Name of phase for refinements to be applied to, or
1186            'all'
1187        """
1188        for proj in self.iter_refinements(refinements, histogram, phase,
1189                                          outputnames, makeBack):
1190            pass
1191        return self
1192
1193    def iter_refinements(self, refinements, histogram='all', phase='all',
1194                         outputnames=None, makeBack=False):
1195        """Conducts a series of refinements, iteratively. Stops after every
1196        refinement and yields this project, to allow error checking or
1197        logging of intermediate results.
1198
1199        >>> def checked_refinements(proj):
1200        ...     for p in proj.iter_refinements(refs):
1201        ...         # Track intermediate results
1202        ...         log(p.histogram('0').residuals)
1203        ...         log(p.phase('0').get_cell())
1204        ...         # Check if parameter diverged, nonsense answer, or whatever
1205        ...         if is_something_wrong(p):
1206        ...             raise Exception("I need a human!")
1207
1208        :param list refinements: A list of dictionaries defining refinements
1209        :param str histogram: Name of histogram for refinements to be applied
1210            to, or 'all'
1211        :param str phase: Name of phase for refinements to be applied to, or
1212            'all'
1213        """
1214        if outputnames:
1215            if len(refinements) != len(outputnames):
1216                raise ValueError("Should have same number of outuputs to"
1217                                 "refinements")
1218        else:
1219            outputnames = [None for r in refinements]
1220
1221        for output, refinement in zip(outputnames, refinements):
1222            self.set_refinement(refinement, histogram)
1223            # Handle 'once' args - refinements that are disabled after this
1224            # refinement
1225            if 'once' in refinement:
1226                temp = {'set': refinement['once']}
1227                self.set_refinement(temp, histogram, phase)
1228
1229            if output:
1230                self.save(output)
1231
1232            self.refine(makeBack=makeBack)
1233            yield self
1234
1235            # Handle 'once' args - refinements that are disabled after this
1236            # refinement
1237            if 'once' in refinement:
1238                temp = {'clear': refinement['once']}
1239                self.set_refinement(temp, histogram, phase)
1240
1241    def set_refinement(self, refinement, histogram='all', phase='all'):
1242        """Apply specified refinements to a given histogram(s) or phase(s).
1243
1244        Refinement parameters are categorize in three groups:
1245
1246        1. Histogram parameters
1247        2. Phase parameters
1248        3. Histogram-and-Phase (HAP) parameters
1249
1250        :param dict refinement: The refinements to be conducted
1251        :param histogram: Either a name of a histogram (str), a list of
1252            histogram names, or 'all' (default)
1253        :param phase: Either a name of a phase (str), a list of phase names, or
1254            'all' (default)
1255
1256        .. seealso::
1257            :func:`~G2PwdrData.set_refinements`
1258            :func:`~G2PwdrData.clear_refinements`
1259            :func:`~G2Phase.set_refinements`
1260            :func:`~G2Phase.clear_refinements`
1261            :func:`~G2Phase.set_HAP_refinements`
1262            :func:`~G2Phase.clear_HAP_refinements`"""
1263
1264        if histogram == 'all':
1265            hists = self.histograms()
1266        elif isinstance(histogram, str) or isinstance(histogram, int):
1267            hists = [self.histogram(histogram)]
1268        else:
1269            hists = [self.histogram(name) for name in histogram]
1270
1271        if phase == 'all':
1272            phases = self.phases()
1273        elif isinstance(phase, str) or isinstance(phase, int):
1274            phases = [self.phase(phase)]
1275        else:
1276            phases = [self.phase(name) for name in phase]
1277
1278        # TODO: HAP parameters:
1279        #   Babinet
1280        #   Extinction
1281        #   HStrain
1282        #   Mustrain
1283        #   Pref. Ori
1284        #   Size
1285
1286        pwdr_set = {}
1287        phase_set = {}
1288        hap_set = {}
1289        for key, val in refinement.get('set', {}).items():
1290            # Apply refinement options
1291            if G2PwdrData.is_valid_refinement_key(key):
1292                pwdr_set[key] = val
1293            elif G2Phase.is_valid_refinement_key(key):
1294                phase_set[key] = val
1295            elif G2Phase.is_valid_HAP_refinement_key(key):
1296                hap_set[key] = val
1297            else:
1298                raise ValueError("Unknown refinement key", key)
1299
1300        for hist in hists:
1301            hist.set_refinements(pwdr_set)
1302        for phase in phases:
1303            phase.set_refinements(phase_set)
1304        for phase in phases:
1305            phase.set_HAP_refinements(hap_set, hists)
1306
1307        pwdr_clear = {}
1308        phase_clear = {}
1309        hap_clear = {}
1310        for key, val in refinement.get('clear', {}).items():
1311            # Clear refinement options
1312            if G2PwdrData.is_valid_refinement_key(key):
1313                pwdr_clear[key] = val
1314            elif G2Phase.is_valid_refinement_key(key):
1315                phase_clear[key] = val
1316            elif G2Phase.is_valid_HAP_refinement_key(key):
1317                hap_set[key] = val
1318            else:
1319                raise ValueError("Unknown refinement key", key)
1320
1321        for hist in hists:
1322            hist.clear_refinements(pwdr_clear)
1323        for phase in phases:
1324            phase.clear_refinements(phase_clear)
1325        for phase in phases:
1326            phase.clear_HAP_refinements(hap_clear, hists)
1327
1328    def index_ids(self):
1329        import GSASIIstrIO as G2strIO
1330        self.save()
1331        return G2strIO.GetUsedHistogramsAndPhases(self.filename)
1332
1333    def add_constraint_raw(self, cons_scope, constr):
1334        """Adds a constraint of type consType to the project.
1335        cons_scope should be one of "Hist", "Phase", "HAP", or "Global".
1336
1337        WARNING it does not check the constraint is well-constructed"""
1338        constrs = self.data['Constraints']['data']
1339        if 'Global' not in constrs:
1340            constrs['Global'] = []
1341        constrs[cons_scope].append(constr)
1342
1343    def hold_many(self, vars, type):
1344        """Apply holds for all the variables in vars, for constraint of a given type.
1345
1346        type is passed directly to add_constraint_raw as consType
1347
1348        :param list vars: A list of variables to hold. Either :class:`GSASIIobj.G2VarObj` objects,
1349            string variable specifiers, or arguments for :meth:`make_var_obj`
1350        :param str type: A string constraint type specifier. See
1351            :class:`~GSASIIscriptable.G2Project.add_constraint_raw`
1352
1353        """
1354        for var in vars:
1355            if isinstance(var, str):
1356                var = self.make_var_obj(var)
1357            elif not isinstance(var, G2obj.G2VarObj):
1358                var = self.make_var_obj(*var)
1359            self.add_constraint_raw(type, [[1.0, var], None, None, 'h'])
1360
1361    def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None,
1362                     reloadIdx=True):
1363        """Wrapper to create a G2VarObj. Takes either a string representaiton ("p:h:name:a")
1364        or individual names of phase, histogram, varname, and atomId.
1365
1366        Automatically converts string phase, hist, or atom names into the ID required
1367        by G2VarObj."""
1368
1369        if reloadIdx:
1370            self.index_ids()
1371
1372        # If string representation, short circuit
1373        if hist is None and varname is None and atomId is None:
1374            if isinstance(phase, str) and ':' in phase:
1375                return G2obj.G2VarObj(phase)
1376
1377        # Get phase index
1378        phaseObj = None
1379        if isinstance(phase, G2Phase):
1380            phaseObj = phase
1381            phase = G2obj.PhaseRanIdLookup[phase.ranId]
1382        elif phase in self.data['Phases']:
1383            phaseObj = self.phase(phase)
1384            phaseRanId = phaseObj.ranId
1385            phase = G2obj.PhaseRanIdLookup[phaseRanId]
1386
1387        # Get histogram index
1388        if isinstance(hist, G2PwdrData):
1389            hist = G2obj.HistRanIdLookup[hist.ranId]
1390        elif hist in self.data:
1391            histRanId = self.histogram(hist).ranId
1392            hist = G2obj.HistRanIdLookup[histRanId]
1393
1394        # Get atom index (if any)
1395        if isinstance(atomId, G2AtomRecord):
1396            atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId]
1397        elif phaseObj:
1398            atomObj = phaseObj.atom(atomId)
1399            if atomObj:
1400                atomRanId = atomObj.ranId
1401                atomId = G2obj.AtomRanIdLookup[phase][atomRanId]
1402
1403        return G2obj.G2VarObj(phase, hist, varname, atomId)
1404
1405
1406class G2AtomRecord(G2ObjectWrapper):
1407    """Wrapper for an atom record. Has convenient accessors via @property.
1408
1409
1410    Available accessors: label, type, refinement_flags, coordinates,
1411        occupancy, ranId, id, adp_flag, uiso
1412
1413    Example:
1414
1415    >>> atom = some_phase.atom("O3")
1416    >>> # We can access the underlying data format
1417    >>> atom.data
1418    ['O3', 'O-2', '', ... ]
1419    >>> # We can also use wrapper accessors
1420    >>> atom.coordinates
1421    (0.33, 0.15, 0.5)
1422    >>> atom.refinement_flags
1423    u'FX'
1424    >>> atom.ranId
1425    4615973324315876477
1426    >>> atom.occupancy
1427    1.0
1428    """
1429    def __init__(self, data, indices, proj):
1430        self.data = data
1431        self.cx, self.ct, self.cs, self.cia = indices
1432        self.proj = proj
1433
1434    @property
1435    def label(self):
1436        return self.data[self.ct-1]
1437
1438    @property
1439    def type(self):
1440        return self.data[self.ct]
1441
1442    @property
1443    def refinement_flags(self):
1444        return self.data[self.ct+1]
1445
1446    @refinement_flags.setter
1447    def refinement_flags(self, other):
1448        # Automatically check it is a valid refinement
1449        for c in other:
1450            if c not in ' FXU':
1451                raise ValueError("Invalid atom refinement: ", other)
1452        self.data[self.ct+1] = unicode(other)
1453
1454    @property
1455    def coordinates(self):
1456        return tuple(self.data[self.cx:self.cx+3])
1457
1458    @property
1459    def occupancy(self):
1460        return self.data[self.cx+3]
1461
1462    @occupancy.setter
1463    def occupancy(self, val):
1464        self.data[self.cx+3] = float(val)
1465
1466    @property
1467    def ranId(self):
1468        return self.data[self.cia+8]
1469
1470    @property
1471    def adp_flag(self):
1472        # Either 'I' or 'A'
1473        return self.data[self.cia]
1474
1475    @property
1476    def uiso(self):
1477        if self.adp_flag == 'I':
1478            return self.data[self.cia+1]
1479        else:
1480            return self.data[self.cia+2:self.cia+8]
1481
1482    @uiso.setter
1483    def uiso(self, value):
1484        if self.adp_flag == 'I':
1485            self.data[self.cia+1] = float(value)
1486        else:
1487            assert len(value) == 6
1488            self.data[self.cia+2:self.cia+8] = [float(v) for v in value]
1489
1490
1491class G2PwdrData(G2ObjectWrapper):
1492    """Wraps a Powder Data Histogram."""
1493    def __init__(self, data, proj):
1494        self.data = data
1495        self.proj = proj
1496
1497    @staticmethod
1498    def is_valid_refinement_key(key):
1499        valid_keys = ['Limits', 'Sample Parameters', 'Background',
1500                      'Instrument Parameters']
1501        return key in valid_keys
1502
1503    @property
1504    def name(self):
1505        return self.data['data'][-1]
1506
1507    @property
1508    def ranId(self):
1509        return self.data['data'][0]['ranId']
1510
1511    @property
1512    def residuals(self):
1513        data = self.data['data'][0]
1514        return {key: data[key]
1515                for key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']}
1516
1517    @property
1518    def id(self):
1519        self.proj.update_ids()
1520        return self.data['data'][0]['hId']
1521
1522    @id.setter
1523    def id(self, val):
1524        self.data['data'][0]['hId'] = val
1525
1526    def fit_fixed_points(self):
1527        """Attempts to apply a background fit to the fixed points currently specified."""
1528        def SetInstParms(Inst):
1529            dataType = Inst['Type'][0]
1530            insVary = []
1531            insNames = []
1532            insVals = []
1533            for parm in Inst:
1534                insNames.append(parm)
1535                insVals.append(Inst[parm][1])
1536                if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1537                    'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1538                        Inst[parm][2] = False
1539            instDict = dict(zip(insNames, insVals))
1540            instDict['X'] = max(instDict['X'], 0.01)
1541            instDict['Y'] = max(instDict['Y'], 0.01)
1542            if 'SH/L' in instDict:
1543                instDict['SH/L'] = max(instDict['SH/L'], 0.002)
1544            return dataType, instDict, insVary
1545
1546        bgrnd = self.data['Background']
1547
1548        # Need our fixed points in order
1549        bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
1550        X = [x for x, y in bgrnd[1]['FixedPoints']]
1551        Y = [y for x, y in bgrnd[1]['FixedPoints']]
1552
1553        limits = self.data['Limits'][1]
1554        if X[0] > limits[0]:
1555            X = [limits[0]] + X
1556            Y = [Y[0]] + Y
1557        if X[-1] < limits[1]:
1558            X += [limits[1]]
1559            Y += [Y[-1]]
1560
1561        # Some simple lookups
1562        controls = self.proj['Controls']['data']
1563        inst, inst2 = self.data['Instrument Parameters']
1564        pwddata = self.data['data'][1]
1565
1566        # Construct the data for background fitting
1567        xBeg = np.searchsorted(pwddata[0], limits[0])
1568        xFin = np.searchsorted(pwddata[0], limits[1])
1569        xdata = pwddata[0][xBeg:xFin]
1570        ydata = si.interp1d(X,Y)(ma.getdata(xdata))
1571
1572        W = [1]*len(xdata)
1573        Z = [0]*len(xdata)
1574
1575        dataType, insDict, insVary = SetInstParms(inst)
1576        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1577
1578        # Do the fit
1579        data = np.array([xdata, ydata, W, Z, Z, Z])
1580        G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
1581                        prevVaryList=bakVary, controls=controls)
1582
1583        # Post-fit
1584        parmDict = {}
1585        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1586        parmDict.update(bakDict)
1587        parmDict.update(insDict)
1588        pwddata[3][xBeg:xFin] *= 0
1589        pwddata[5][xBeg:xFin] *= 0
1590        pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
1591
1592        # TODO adjust pwddata? GSASIIpwdGUI.py:1041
1593        # TODO update background
1594        self.proj.save()
1595
1596    def y_calc(self):
1597        return self.data['data'][1][3]
1598
1599    def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True):
1600        try:
1601            import matplotlib.pyplot as plt
1602            data = self.data['data'][1]
1603            if Yobs:
1604                plt.plot(data[0], data[1], label='Yobs')
1605            if Ycalc:
1606                plt.plot(data[0], data[3], label='Ycalc')
1607            if Background:
1608                plt.plot(data[0], data[4], label='Background')
1609            if Residual:
1610                plt.plot(data[0], data[5], label="Residual")
1611        except ImportError:
1612            pass
1613
1614    def set_refinements(self, refs):
1615        """Sets the refinement parameter 'key' to the specification 'value'
1616
1617        :param dict refs: A dictionary of the parameters to be set. See
1618                          :ref:`Histogram_parameters_table` for a description of
1619                          what these dictionaries should be.
1620
1621        :returns: None
1622
1623        """
1624        do_fit_fixed_points = False
1625        for key, value in refs.items():
1626            if key == 'Limits':
1627                old_limits = self.data['Limits'][1]
1628                new_limits = value
1629                if isinstance(new_limits, dict):
1630                    if 'low' in new_limits:
1631                        old_limits[0] = new_limits['low']
1632                    if 'high' in new_limits:
1633                        old_limits[1] = new_limits['high']
1634                else:
1635                    old_limits[0], old_limits[1] = new_limits
1636            elif key == 'Sample Parameters':
1637                sample = self.data['Sample Parameters']
1638                for sparam in value:
1639                    if sparam not in sample:
1640                        raise ValueError("Unknown refinement parameter, "
1641                                         + str(sparam))
1642                    sample[sparam][1] = True
1643            elif key == 'Background':
1644                bkg, peaks = self.data['Background']
1645
1646                # If True or False, just set the refine parameter
1647                if value in (True, False):
1648                    bkg[1] = value
1649                    return
1650
1651                if 'type' in value:
1652                    bkg[0] = value['type']
1653                if 'refine' in value:
1654                    bkg[1] = value['refine']
1655                if 'no. coeffs' in value:
1656                    cur_coeffs = bkg[2]
1657                    n_coeffs = value['no. coeffs']
1658                    if n_coeffs > cur_coeffs:
1659                        for x in range(n_coeffs - cur_coeffs):
1660                            bkg.append(0.0)
1661                    else:
1662                        for _ in range(cur_coeffs - n_coeffs):
1663                            bkg.pop()
1664                    bkg[2] = n_coeffs
1665                if 'coeffs' in value:
1666                    bkg[3:] = value['coeffs']
1667                if 'FixedPoints' in value:
1668                    peaks['FixedPoints'] = [(float(a), float(b))
1669                                            for a, b in value['FixedPoints']]
1670                if value.get('fit fixed points', False):
1671                    do_fit_fixed_points = True
1672
1673            elif key == 'Instrument Parameters':
1674                instrument, secondary = self.data['Instrument Parameters']
1675                for iparam in value:
1676                    try:
1677                        instrument[iparam][2] = True
1678                    except IndexError:
1679                        raise ValueError("Invalid key:", iparam)
1680            else:
1681                raise ValueError("Unknown key:", key)
1682        # Fit fixed points after the fact - ensure they are after fixed points
1683        # are added
1684        if do_fit_fixed_points:
1685            # Background won't be fit if refinement flag not set
1686            orig = self.data['Background'][0][1]
1687            self.data['Background'][0][1] = True
1688            self.fit_fixed_points()
1689            # Restore the previous value
1690            self.data['Background'][0][1] = orig
1691
1692    def clear_refinements(self, refs):
1693        """Clears the refinement parameter 'key' and its associated value.
1694
1695        :param dict refs: A dictionary of parameters to clear."""
1696        for key, value in refs.items():
1697            if key == 'Limits':
1698                old_limits, cur_limits = self.data['Limits']
1699                cur_limits[0], cur_limits[1] = old_limits
1700            elif key == 'Sample Parameters':
1701                sample = self.data['Sample Parameters']
1702                for sparam in value:
1703                    sample[sparam][1] = False
1704            elif key == 'Background':
1705                bkg, peaks = self.data['Background']
1706
1707                # If True or False, just set the refine parameter
1708                if value in (True, False):
1709                    bkg[1] = False
1710                    return
1711
1712                bkg[1] = False
1713                if 'FixedPoints' in value:
1714                    if 'FixedPoints' in peaks:
1715                        del peaks['FixedPoints']
1716            elif key == 'Instrument Parameters':
1717                instrument, secondary = self.data['Instrument Parameters']
1718                for iparam in value:
1719                    instrument[iparam][2] = False
1720            else:
1721                raise ValueError("Unknown key:", key)
1722
1723
1724class G2Phase(G2ObjectWrapper):
1725    """A wrapper object around a given phase.
1726
1727    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1728    """
1729    def __init__(self, data, name, proj):
1730        self.data = data
1731        self.name = name
1732        self.proj = proj
1733
1734    @staticmethod
1735    def is_valid_refinement_key(key):
1736        valid_keys = ["Cell", "Atoms", "LeBail"]
1737        return key in valid_keys
1738
1739    @staticmethod
1740    def is_valid_HAP_refinement_key(key):
1741        valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain",
1742                      "Pref.Ori.", "Show", "Size", "Use", "Scale"]
1743        return key in valid_keys
1744
1745    def atom(self, atomlabel):
1746        """Returns the atom specified by atomlabel, or None if it does not
1747        exist.
1748
1749        :param str atomlabel: The name of the atom (e.g. "O2")
1750        :returns: A :class:`G2AtomRecord` object
1751            representing the atom.
1752        """
1753        # Consult GSASIIobj.py for the meaning of this
1754        cx, ct, cs, cia = self.data['General']['AtomPtrs']
1755        ptrs = [cx, ct, cs, cia]
1756        atoms = self.data['Atoms']
1757        for atom in atoms:
1758            if atom[ct-1] == atomlabel:
1759                return G2AtomRecord(atom, ptrs, self.proj)
1760
1761    def atoms(self):
1762        """Returns a list of atoms present in the phase.
1763
1764        :returns: A list of :class:`G2AtomRecord` objects.
1765
1766        .. seealso::
1767            :func:`~GSASIIscriptable.G2Phase.atom`
1768            :class:`G2AtomRecord`
1769        """
1770        ptrs = self.data['General']['AtomPtrs']
1771        output = []
1772        atoms = self.data['Atoms']
1773        for atom in atoms:
1774            output.append(G2AtomRecord(atom, ptrs, self.proj))
1775        return output
1776
1777    def histograms(self):
1778        output = []
1779        for hname in self.data.get('Histograms', {}).keys():
1780            output.append(self.proj.histogram(hname))
1781        return output
1782
1783    @property
1784    def ranId(self):
1785        return self.data['ranId']
1786
1787    @property
1788    def id(self):
1789        return self.data['pId']
1790
1791    @id.setter
1792    def id(self, val):
1793        self.data['pId'] = val
1794
1795    def get_cell(self):
1796        """Returns a dictionary of the cell parameters, with keys:
1797            'length_a', 'length_b', 'length_c', 'angle_alpha', 'angle_beta', 'angle_gamma', 'volume'
1798
1799        :returns: a dict
1800
1801        .. seealso::
1802           :func:`~G2Phase.get_cell_and_esd`
1803
1804        """
1805        cell = self.data['General']['Cell']
1806        return {'length_a': cell[1], 'length_b': cell[2], 'length_c': cell[3],
1807                'angle_alpha': cell[4], 'angle_beta': cell[5], 'angle_gamma': cell[6],
1808                'volume': cell[7]}
1809
1810    def get_cell_and_esd(self):
1811        """
1812        Returns a pair of dictionaries, the first representing the unit cell, the second
1813        representing the estimated standard deviations of the unit cell.
1814
1815        :returns: a tuple of two dictionaries
1816
1817        .. seealso::
1818           :func:`~G2Phase.get_cell`
1819
1820        """
1821        # translated from GSASIIstrIO.ExportBaseclass.GetCell
1822        import GSASIIstrIO as G2stIO
1823        import GSASIIlattice as G2lat
1824        import GSASIImapvars as G2mv
1825        try:
1826            pfx = str(self.id) + '::'
1827            sgdata = self['General']['SGData']
1828            covDict = self.proj['Covariance']['data']
1829
1830            parmDict = dict(zip(covDict.get('varyList',[]),
1831                                covDict.get('variables',[])))
1832            sigDict = dict(zip(covDict.get('varyList',[]),
1833                               covDict.get('sig',[])))
1834
1835            if covDict.get('covMatrix') is not None:
1836                sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
1837                                                  covDict['varyList'],
1838                                                  parmDict))
1839
1840            A, sigA = G2stIO.cellFill(pfx, sgdata, parmDict, sigDict)
1841            cellSig = G2stIO.getCellEsd(pfx, sgdata, A, self.proj['Covariance']['data'])
1842            cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),)
1843            cellDict, cellSigDict = {}, {}
1844            for i, key in enumerate(['length_a', 'length_b', 'length_c',
1845                                     'angle_alpha', 'angle_beta', 'angle_gamma',
1846                                     'volume']):
1847                cellDict[key] = cellList[i]
1848                cellSigDict[key] = cellSig[i]
1849            return cellDict, cellSigDict
1850        except KeyError:
1851            cell = self.get_cell()
1852            return cell, {key: 0.0 for key in cell}
1853
1854    def export_CIF(self, outputname, quickmode=True):
1855        """Write this phase to a .cif file named outputname
1856
1857        :param str outputname: The name of the .cif file to write to
1858        :param bool quickmode: Currently ignored. Carryover from exports.G2export_CIF"""
1859        # This code is all taken from exports/G2export_CIF.py
1860        # Functions copied have the same names
1861        import GSASIImath as G2mth
1862        import GSASIImapvars as G2mv
1863        from exports import G2export_CIF as cif
1864
1865        CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M")
1866        CIFname = os.path.splitext(self.proj.filename)[0]
1867        CIFname = os.path.split(CIFname)[1]
1868        CIFname = ''.join([c if ord(c) < 128 else ''
1869                           for c in CIFname.replace(' ', '_')])
1870        try:
1871            author = self.proj['Controls']['data'].get('Author','').strip()
1872        except KeyError:
1873            pass
1874        oneblock = True
1875
1876        covDict = self.proj['Covariance']['data']
1877        parmDict = dict(zip(covDict.get('varyList',[]),
1878                            covDict.get('variables',[])))
1879        sigDict = dict(zip(covDict.get('varyList',[]),
1880                           covDict.get('sig',[])))
1881
1882        if covDict.get('covMatrix') is not None:
1883            sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
1884                                              covDict['varyList'],
1885                                              parmDict))
1886
1887        with open(outputname, 'w') as fp:
1888            fp.write(' \n' + 70*'#' + '\n')
1889            cif.WriteCIFitem(fp, 'data_' + CIFname)
1890            # from exports.G2export_CIF.WritePhaseInfo
1891            cif.WriteCIFitem(fp, '\n# phase info for '+str(self.name) + ' follows')
1892            cif.WriteCIFitem(fp, '_pd_phase_name', self.name)
1893            # TODO get esds
1894            cellDict = self.get_cell()
1895            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
1896            names = ['length_a','length_b','length_c',
1897                     'angle_alpha','angle_beta ','angle_gamma',
1898                     'volume']
1899            for key, val in cellDict.items():
1900                cif.WriteCIFitem(fp, '_cell_' + key, G2mth.ValEsd(val))
1901
1902            cif.WriteCIFitem(fp, '_symmetry_cell_setting',
1903                         self.data['General']['SGData']['SGSys'])
1904
1905            spacegroup = self.data['General']['SGData']['SpGrp'].strip()
1906            # regularize capitalization and remove trailing H/R
1907            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
1908            cif.WriteCIFitem(fp, '_symmetry_space_group_name_H-M', spacegroup)
1909
1910            # generate symmetry operations including centering and center of symmetry
1911            SymOpList, offsetList, symOpList, G2oprList, G2opcodes = G2spc.AllOps(
1912                self.data['General']['SGData'])
1913            cif.WriteCIFitem(fp, 'loop_\n    _space_group_symop_id\n    _space_group_symop_operation_xyz')
1914            for i, op in enumerate(SymOpList,start=1):
1915                cif.WriteCIFitem(fp, '   {:3d}  {:}'.format(i,op.lower()))
1916
1917            # TODO skipped histograms, exports/G2export_CIF.py:880
1918
1919            # report atom params
1920            if self.data['General']['Type'] in ['nuclear','macromolecular']:        #this needs macromolecular variant, etc!
1921                cif.WriteAtomsNuclear(fp, self.data, self.name, parmDict, sigDict, [])
1922                # self._WriteAtomsNuclear(fp, parmDict, sigDict)
1923            else:
1924                raise Exception("no export for "+str(self.data['General']['Type'])+" coordinates implemented")
1925            # report cell contents
1926            cif.WriteComposition(fp, self.data, self.name, parmDict)
1927            if not quickmode and self.data['General']['Type'] == 'nuclear':      # report distances and angles
1928                # WriteDistances(fp,self.name,SymOpList,offsetList,symOpList,G2oprList)
1929                raise NotImplementedError("only quickmode currently supported")
1930            if 'Map' in self.data['General'] and 'minmax' in self.data['General']['Map']:
1931                cif.WriteCIFitem(fp,'\n# Difference density results')
1932                MinMax = self.data['General']['Map']['minmax']
1933                cif.WriteCIFitem(fp,'_refine_diff_density_max',G2mth.ValEsd(MinMax[0],-0.009))
1934                cif.WriteCIFitem(fp,'_refine_diff_density_min',G2mth.ValEsd(MinMax[1],-0.009))
1935
1936
1937    def set_refinements(self, refs):
1938        """Sets the refinement parameter 'key' to the specification 'value'
1939
1940        :param dict refs: A dictionary of the parameters to be set. See
1941                          :ref:`Phase_parameters_table` for a description of
1942                          this dictionary.
1943
1944        :returns: None"""
1945        for key, value in refs.items():
1946            if key == "Cell":
1947                self.data['General']['Cell'][0] = value
1948
1949            elif key == "Atoms":
1950                for atomlabel, atomrefinement in value.items():
1951                    if atomlabel == 'all':
1952                        for atom in self.atoms():
1953                            atom.refinement_flags = atomrefinement
1954                    else:
1955                        atom = self.atom(atomlabel)
1956                        if atom is None:
1957                            raise ValueError("No such atom: " + atomlabel)
1958                        atom.refinement_flags = atomrefinement
1959
1960            elif key == "LeBail":
1961                hists = self.data['Histograms']
1962                for hname, hoptions in hists.items():
1963                    if 'LeBail' not in hoptions:
1964                        hoptions['newLeBail'] = bool(True)
1965                    hoptions['LeBail'] = bool(value)
1966            else:
1967                raise ValueError("Unknown key:", key)
1968
1969    def clear_refinements(self, refs):
1970        """Clears a given set of parameters.
1971
1972        :param dict refs: The parameters to clear"""
1973        for key, value in refs.items():
1974            if key == "Cell":
1975                self.data['General']['Cell'][0] = False
1976            elif key == "Atoms":
1977                cx, ct, cs, cia = self.data['General']['AtomPtrs']
1978
1979                for atomlabel in value:
1980                    atom = self.atom(atomlabel)
1981                    # Set refinement to none
1982                    atom.refinement_flags = ' '
1983            elif key == "LeBail":
1984                hists = self.data['Histograms']
1985                for hname, hoptions in hists.items():
1986                    if 'LeBail' not in hoptions:
1987                        hoptions['newLeBail'] = True
1988                    hoptions['LeBail'] = False
1989            else:
1990                raise ValueError("Unknown key:", key)
1991
1992    def set_HAP_refinements(self, refs, histograms='all'):
1993        """Sets the given HAP refinement parameters between this phase and
1994        the given histograms
1995
1996        :param dict refs: A dictionary of the parameters to be set. See
1997                          :ref:`HAP_parameters_table` for a description of this
1998                          dictionary.
1999        :param histograms: Either 'all' (default) or a list of the histograms
2000            whose HAP parameters will be set with this phase. Histogram and phase
2001            must already be associated
2002
2003        :returns: None
2004        """
2005        if histograms == 'all':
2006            histograms = self.data['Histograms'].values()
2007        else:
2008            histograms = [self.data['Histograms'][h.name] for h in histograms]
2009
2010        for key, val in refs.items():
2011            for h in histograms:
2012                if key == 'Babinet':
2013                    try:
2014                        sets = list(val)
2015                    except ValueError:
2016                        sets = ['BabA', 'BabU']
2017                    for param in sets:
2018                        if param not in ['BabA', 'BabU']:
2019                            raise ValueError("Not sure what to do with" + param)
2020                        for hist in histograms:
2021                            hist['Babinet'][param][1] = True
2022                elif key == 'Extinction':
2023                    for h in histograms:
2024                        h['Extinction'][1] = bool(val)
2025                elif key == 'HStrain':
2026                    for h in histograms:
2027                        hist['HStrain'][1] = [bool(val) for p in hist['Hstrain'][0]]
2028                elif key == 'Mustrain':
2029                    for h in histograms:
2030                        mustrain = h['Mustrain']
2031                        newType = None
2032                        if isinstance(val, (unicode, str)):
2033                            if val in ['isotropic', 'uniaxial', 'generalized']:
2034                                newType = val
2035                            else:
2036                                raise ValueError("Not a Mustrain type: " + val)
2037                        elif isinstance(val, dict):
2038                            newType = val.get('type', None)
2039                            direction = val.get('direction', None)
2040
2041                        if newType:
2042                            mustrain[0] = newType
2043                            if newType == 'isotropic':
2044                                mustrain[2][0] = True
2045                                mustrain[5] = [False for p in mustrain[4]]
2046                            elif newType == 'uniaxial':
2047                                if 'refine' in val:
2048                                    types = val['refine']
2049                                    if isinstance(types, (unicode, str)):
2050                                        types = [types]
2051                                    elif isinstance(types, bool):
2052                                        mustrain[2][0] = types
2053                                        mustrain[2][1] = types
2054                                        types = []
2055                                    else:
2056                                        raise ValueError("Not sure what to do with: "
2057                                                         + str(types))
2058                                else:
2059                                    types = []
2060
2061                                for unitype in types:
2062                                    if unitype == 'equatorial':
2063                                        mustrain[2][0] = True
2064                                    elif unitype == 'axial':
2065                                        mustrain[2][1] = True
2066                                    else:
2067                                        msg = 'Invalid uniaxial mustrain type'
2068                                        raise ValueError(msg + ': ' + unitype)
2069                            else:  # newtype == 'generalized'
2070                                mustrain[2] = [False for p in mustrain[1]]
2071
2072                        if direction:
2073                            direction = [int(n) for n in direction]
2074                            if len(direction) != 3:
2075                                raise ValueError("Expected hkl, found", direction)
2076                            mustrain[3] = direction
2077                elif key == 'Pref.Ori.':
2078                    for h in histograms:
2079                        h['Pref.Ori.'][2] = bool(val)
2080                elif key == 'Show':
2081                    for h in histograms:
2082                        h['Show'] = bool(val)
2083                elif key == 'Size':
2084                    # TODO
2085                    raise NotImplementedError()
2086                elif key == 'Use':
2087                    for h in histograms:
2088                        h['Use'] = bool(val)
2089                elif key == 'Scale':
2090                    for h in histograms:
2091                        h['Scale'][1] = bool(val)
2092
2093    def clear_HAP_refinements(self, refs, histograms='all'):
2094        """Clears the given HAP refinement parameters between this phase and
2095        the given histograms
2096
2097        :param dict refs: A dictionary of the parameters to be cleared.
2098        :param histograms: Either 'all' (default) or a list of the histograms
2099            whose HAP parameters will be cleared with this phase. Histogram and
2100            phase must already be associated
2101
2102        :returns: None
2103        """
2104        if histograms == 'all':
2105            histograms = self.data['Histograms'].values()
2106        else:
2107            histograms = [self.data['Histograms'][h.name] for h in histograms]
2108
2109        for key, val in refs.items():
2110            for h in histograms:
2111                if key == 'Babinet':
2112                    try:
2113                        sets = list(val)
2114                    except ValueError:
2115                        sets = ['BabA', 'BabU']
2116                    for param in sets:
2117                        if param not in ['BabA', 'BabU']:
2118                            raise ValueError("Not sure what to do with" + param)
2119                        for hist in histograms:
2120                            hist['Babinet'][param][1] = False
2121                elif key == 'Extinction':
2122                    for h in histograms:
2123                        h['Extinction'][1] = False
2124                elif key == 'HStrain':
2125                    for h in histograms:
2126                        hist['HStrain'][1] = [False for p in hist['Hstrain'][0]]
2127                elif key == 'Mustrain':
2128                    for h in histograms:
2129                        mustrain = h['Mustrain']
2130                        mustrain[2] = [False for p in mustrain[1]]
2131                        mustrain[5] = [False for p in mustrain[4]]
2132                elif key == 'Pref.Ori.':
2133                    for h in histograms:
2134                        h['Pref.Ori.'][2] = False
2135                elif key == 'Show':
2136                    for h in histograms:
2137                        h['Show'] = False
2138                elif key == 'Size':
2139                    raise NotImplementedError()
2140                elif key == 'Use':
2141                    for h in histograms:
2142                        h['Use'] = False
2143                elif key == 'Scale':
2144                    for h in histograms:
2145                        h['Scale'][1] = False
2146
2147
2148##########################
2149# Command Line Interface #
2150##########################
2151# Each of these takes an argparse.Namespace object as their argument,
2152# representing the parsed command-line arguments for the relevant subcommand.
2153# The argument specification for each is in the subcommands dictionary (see
2154# below)
2155
2156
2157def create(args):
2158    """The create subcommand."""
2159    proj = G2Project(filename=args.filename)
2160
2161    hist_objs = []
2162    for h in args.histograms:
2163        hist_objs.append(proj.add_powder_histogram(h, args.iparams))
2164
2165    for p in args.phases:
2166        proj.add_phase(p, histograms=hist_objs)
2167
2168    proj.save()
2169
2170
2171def dump(args):
2172    """The dump subcommand"""
2173    if not args.histograms and not args.phases:
2174        args.raw = True
2175    if args.raw:
2176        import IPython.lib.pretty as pretty
2177
2178    for fname in args.files:
2179        if args.raw:
2180            proj, nameList = LoadDictFromProjFile(fname)
2181            print("file:", fname)
2182            print("names:", nameList)
2183            for key, val in proj.items():
2184                print(key, ":")
2185                pretty.pprint(val)
2186        else:
2187            proj = G2Project(fname)
2188            if args.histograms:
2189                hists = proj.histograms()
2190                for h in hists:
2191                    print(fname, "hist", h.id, h.name)
2192            if args.phases:
2193                phase_list = proj.phases()
2194                for p in phase_list:
2195                    print(fname, "phase", p.id, p.name)
2196
2197
2198def IPyBrowse(args):
2199    """Load a .gpx file and then open a IPython shell to browse it
2200    """
2201    for fname in args.files:
2202        proj, nameList = LoadDictFromProjFile(fname)
2203        msg = "\nfname {} loaded into proj (dict) with names in nameList".format(fname)
2204        GSASIIpath.IPyBreak_base(msg)
2205        break
2206
2207
2208def refine(args):
2209    """The refine subcommand"""
2210    proj = G2Project(args.gpxfile)
2211    if args.refinements is None:
2212        proj.refine()
2213    else:
2214        import json
2215        with open(args.refinements) as refs:
2216            refs = json.load(refs)
2217        proj.do_refinements(refs['refinements'])
2218
2219
2220def seqrefine(args):
2221    """The seqrefine subcommand"""
2222    raise NotImplementedError("seqrefine is not yet implemented")
2223
2224
2225def export(args):
2226    """The export subcommand"""
2227    # Export phase as CIF to args.exportfile
2228    proj = G2Project(args.gpxfile)
2229    phase = proj.phase(args.phase)
2230    phase.export_CIF(args.exportfile)
2231
2232
2233def _args_kwargs(*args, **kwargs):
2234    return args, kwargs
2235
2236# A dictionary of the name of each subcommand, and a tuple
2237# of its associated function and a list of its arguments
2238# The arguments are passed directly to the add_argument() method
2239# of an argparse.ArgumentParser
2240subcommands = {"create":
2241               (create, [_args_kwargs('filename',
2242                                      help='the project file to create. should end in .gpx'),
2243
2244                         _args_kwargs('-i', '--iparams',
2245                                      help='instrument parameter file'),
2246
2247                         _args_kwargs('-d', '--histograms',
2248                                      nargs='+',
2249                                      help='list of datafiles to add as histograms'),
2250
2251                         _args_kwargs('-p', '--phases',
2252                                      nargs='+',
2253                                      help='list of phases to add. phases are '
2254                                           'automatically associated with all '
2255                                           'histograms given.')]),
2256
2257               "dump": (dump, [_args_kwargs('-d', '--histograms',
2258                                     action='store_true',
2259                                     help='list histograms in files, overrides --raw'),
2260
2261                               _args_kwargs('-p', '--phases',
2262                                            action='store_true',
2263                                            help='list phases in files, overrides --raw'),
2264
2265                               _args_kwargs('-r', '--raw',
2266                                      action='store_true', help='dump raw file contents, default'),
2267
2268                               _args_kwargs('files', nargs='+')]),
2269
2270               "refine":
2271               (refine, [_args_kwargs('gpxfile', help='the project file to refine'),
2272                         _args_kwargs('refinements',
2273                                      help='json file of refinements to apply. if not present'
2274                                           ' refines file as-is',
2275                                      default=None,
2276                                      nargs='?')]),
2277
2278               "seqrefine": (seqrefine, []),
2279               "export": (export, [_args_kwargs('gpxfile',
2280                                                help='the project file from which to export'),
2281                                   _args_kwargs('phase', help='identifier of phase to export'),
2282                                   _args_kwargs('exportfile', help='the .cif file to export to')]),
2283               "browse": (IPyBrowse, [_args_kwargs('files', nargs='+',
2284                                                   help='list of files to browse')])}
2285
2286
2287def main():
2288    '''The command line interface for GSASIIscriptable.
2289
2290    Executes one of the following subcommands:
2291
2292        * :func:`create`
2293        * :func:`dump`
2294        * :func:`refine`
2295        * :func:`seqrefine`
2296        * :func:`export`
2297        * :func:`browse`
2298
2299    .. seealso::
2300        :func:`create`
2301        :func:`dump`
2302        :func:`refine`
2303        :func:`seqrefine`
2304        :func:`export`
2305        :func:`browse`
2306    '''
2307    parser = argparse.ArgumentParser()
2308    subs = parser.add_subparsers()
2309
2310    # Create all of the specified subparsers
2311    for name, (func, args) in subcommands.items():
2312        new_parser = subs.add_parser(name)
2313        for listargs, kwds in args:
2314            new_parser.add_argument(*listargs, **kwds)
2315        new_parser.set_defaults(func=func)
2316
2317    # Parse and trigger subcommand
2318    result = parser.parse_args()
2319    result.func(result)
2320
2321if __name__ == '__main__':
2322    main()
Note: See TracBrowser for help on using the repository browser.