source: trunk/GSASIIscriptable.py @ 3102

Last change on this file since 3102 was 3102, checked in by toby, 8 years ago

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